[GRASS-SVN] r70973 - in grass-addons/grass7/raster: . r.sim.water.mp
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Apr 27 21:43:41 PDT 2017
Author: jhofier
Date: 2017-04-27 21:43:41 -0700 (Thu, 27 Apr 2017)
New Revision: 70973
Added:
grass-addons/grass7/raster/r.sim.water.mp/
grass-addons/grass7/raster/r.sim.water.mp/Makefile
grass-addons/grass7/raster/r.sim.water.mp/description-pak
grass-addons/grass7/raster/r.sim.water.mp/erod.c
grass-addons/grass7/raster/r.sim.water.mp/hydro.c
grass-addons/grass7/raster/r.sim.water.mp/input.c
grass-addons/grass7/raster/r.sim.water.mp/main.c
grass-addons/grass7/raster/r.sim.water.mp/observation_points.c
grass-addons/grass7/raster/r.sim.water.mp/output.c
grass-addons/grass7/raster/r.sim.water.mp/r.sim.water.mp.html
grass-addons/grass7/raster/r.sim.water.mp/r_sim_water.png
grass-addons/grass7/raster/r.sim.water.mp/random.c
grass-addons/grass7/raster/r.sim.water.mp/simlib.h
grass-addons/grass7/raster/r.sim.water.mp/spearfish.sh
grass-addons/grass7/raster/r.sim.water.mp/utils.c
grass-addons/grass7/raster/r.sim.water.mp/waterglobs.h
Log:
r.sim.water with openmp suppport
Added: grass-addons/grass7/raster/r.sim.water.mp/Makefile
===================================================================
--- grass-addons/grass7/raster/r.sim.water.mp/Makefile (rev 0)
+++ grass-addons/grass7/raster/r.sim.water.mp/Makefile 2017-04-28 04:43:41 UTC (rev 70973)
@@ -0,0 +1,17 @@
+
+MODULE_TOPDIR = ../../..
+
+PGM=r.sim.water.mp
+
+EXTRA_CLEAN_DIRS=doxygenhtml
+
+LIBES = $(GISLIB) $(DATETIMELIB) $(GMATHLIB) $(RASTERLIB) $(VECTORLIB)
+DEPENDENCIES = $(GISDEP) $(BITMAPDEP) $(DBMIDEP) $(GMATHDEP) $(LINKMDEP) $(XDRDEP)
+EXTRA_INC = $(VECT_INC)
+EXTRA_CFLAGS= $(VECT_CFLAGS) -fopenmp
+EXTRA_LIBS = $(OCLLIB) $(GISLIB) -lgomp $(MATHLIB)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
+
Added: grass-addons/grass7/raster/r.sim.water.mp/description-pak
===================================================================
--- grass-addons/grass7/raster/r.sim.water.mp/description-pak (rev 0)
+++ grass-addons/grass7/raster/r.sim.water.mp/description-pak 2017-04-28 04:43:41 UTC (rev 70973)
@@ -0,0 +1 @@
+Package created with checkinstall 1.6.2
Added: grass-addons/grass7/raster/r.sim.water.mp/erod.c
===================================================================
--- grass-addons/grass7/raster/r.sim.water.mp/erod.c (rev 0)
+++ grass-addons/grass7/raster/r.sim.water.mp/erod.c 2017-04-28 04:43:41 UTC (rev 70973)
@@ -0,0 +1,68 @@
+
+/****************************************************************************
+ *
+ * MODULE: simwe library
+ * AUTHOR(S): Helena Mitasova, Jaro Hofierka, Lubos Mitas:
+ * PURPOSE: Hydrologic and sediment transport simulation (SIMWE)
+ *
+ * COPYRIGHT: (C) 2002 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.
+ *
+ *****************************************************************************/
+
+/* erod.c (simlib), 20.nov.2002, JH */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <grass/gis.h>
+#include <grass/bitmap.h>
+#include <grass/linkm.h>
+
+#include "waterglobs.h"
+
+/* divergence computation from a given field */
+
+void erod(double **hw)
+{
+ /* hw = sigma or gamma */
+
+ double dyp, dyn, dya, dxp, dxn, dxa;
+ int k, l;
+ int l1, lp, k1, kp, ln, kn, k2, l2;
+
+
+ for (k = 0; k < my; k++) {
+ for (l = 0; l < mx; l++) {
+ lp = max(0, l - 2);
+ l1 = lp + 1;
+ kp = max(0, k - 2);
+ k1 = kp + 1;
+ ln = min(mx - 1, l + 1);
+ l2 = ln - 1;
+ kn = min(my - 1, k + 1);
+ k2 = kn - 1;
+
+ if (zz[k][l] != UNDEF || zz[k][ln] != UNDEF || zz[kp][l] != UNDEF || zz[k][lp] != UNDEF || zz[k][l1] != UNDEF || zz[k1][l] != UNDEF || zz[kn][l] != UNDEF) { /* jh fix */
+
+ dxp = (v1[k][lp] * hw[k][lp] - v1[k][l1] * hw[k][l1]) / stepx;
+ dxn = (v1[k][l2] * hw[k][l2] - v1[k][ln] * hw[k][ln]) / stepx;
+ dxa = 0.5 * (dxp + dxn);
+
+
+ dyp = (v2[kp][l] * hw[kp][l] - v2[k1][l] * hw[k1][l]) / stepy;
+ dyn = (v2[k2][l] * hw[k2][l] - v2[kn][l] * hw[kn][l]) / stepy;
+ dya = 0.5 * (dyp + dyn);
+
+ er[k][l] = (dxa + dya) / deltap;
+
+ }
+ else
+ er[k][l] = UNDEF;
+
+ }
+ }
+}
Added: grass-addons/grass7/raster/r.sim.water.mp/hydro.c
===================================================================
--- grass-addons/grass7/raster/r.sim.water.mp/hydro.c (rev 0)
+++ grass-addons/grass7/raster/r.sim.water.mp/hydro.c 2017-04-28 04:43:41 UTC (rev 70973)
@@ -0,0 +1,619 @@
+
+/****************************************************************************
+ *
+ * MODULE: simwe library
+ * AUTHOR(S): Helena Mitasova, Jaro Hofierka, Lubos Mitas:
+ * PURPOSE: Hydrologic and sediment transport simulation (SIMWE)
+ *
+ * COPYRIGHT: (C) 2002 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.
+ *
+ *****************************************************************************/
+
+/* hydro.c (simlib), 20.nov.2002, JH */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <grass/gis.h>
+#include <grass/bitmap.h>
+#include <grass/linkm.h>
+#include <grass/glocale.h>
+
+#include "waterglobs.h"
+#include "simlib.h"
+
+/*
+ * Soeren 8. Mar 2011 TODO:
+ * Put all these global variables into several meaningful structures and
+ * document use and purpose.
+ *
+ */
+
+struct _points points;
+
+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;
+
+struct seed seed;
+
+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 bxmi, bymi;
+
+float **zz, **cchez;
+double **v1, **v2, **slope;
+double **gama, **gammas, **si, **inf, **sigma;
+float **dc, **tau, **er, **ct, **trap;
+float **dif;
+
+/* suspected BUG below: fist array subscripts go from 1 to MAXW
+ * correct: from 0 to MAXW - 1, e.g for (lw = 0; lw < MAXW; lw++) */
+double vavg[MAXW][2], stack[MAXW][3], w[MAXW][3];
+int iflag[MAXW];
+
+double hbeta;
+double hhmax, sisum, vmean;
+double infsum, infmean;
+int maxw, maxwa, nwalk;
+double rwalk, bresx, bresy, xrand, yrand;
+double stepx, stepy, xp0, yp0;
+double chmean, si0, deltap, deldif, cch, hhc, halpha;
+double eps;
+int maxwab, nstack;
+int iterout, mx2o, my2o;
+int miter, nwalka;
+double timec;
+int ts, timesec;
+
+double rain_val;
+double manin_val;
+double infil_val;
+
+struct History history; /* holds meta-data (title, comments,..) */
+
+/* **************************************************** */
+/* create walker representation of si */
+/* ******************************************************** */
+/* .......... iblock loop */
+
+void main_loop(void)
+{
+ G_message(_("IT'S ALIVE"));
+ G_message(_("Starting MAINLOOP"));
+#ifdef PARALLEL
+ printTimeDiff("MP");
+#endif
+
+ int i, ii, l, k;
+ int icoub, nmult;
+ int iw, iblock, lw;
+ int itime, iter1;
+ int nfiterh, nfiterw;
+ int mgen, mgen2, mgen3;
+ int nblock;
+ int icfl;
+ int mitfac;
+/* int mitfac, p; */
+ double x, y;
+ double velx, vely, stxm, stym;
+ double factor, conn, gaux, gauy;
+ double d1, addac, decr;
+ double barea, sarea, walkwe;
+ double gen, gen2, wei2, wei3, wei, weifac;
+ float eff;
+
+ nblock = 1;
+ icoub = 0;
+ icfl = 0;
+ nstack = 0;
+
+ if (maxwa > (MAXW - mx * my)) {
+ mitfac = maxwa / (MAXW - mx * my);
+ nblock = mitfac + 1;
+ maxwa = maxwa / nblock;
+ }
+
+ /* Create the observation points */
+ create_observation_points();
+
+ G_debug(2, " maxwa, nblock %d %d", maxwa, nblock);
+
+//---------------------------------------------------------------------------------------------------------------
+
+#ifdef PARALLEL
+ G_message(_("Starting parallel work"));
+ printTimeDiff("L0");
+#endif
+
+ int cnt1, cnt2;
+ /* int **lw_start = NULL;
+ if (lw_start == NULL) {
+ lw_start = (int **)G_calloc((my), sizeof(int **));
+ for (cnt1 = 0; cnt1 < my; cnt1++) {
+ lw_start[cnt1] = (int *)G_calloc((mx), sizeof(int *));
+ }
+
+
+ }*/
+
+ for (iblock = 1; iblock <= nblock; iblock++) {
+ ++icoub;
+
+ lw = 0;
+ walkwe = 0.;
+ barea = stepx * stepy;
+ sarea = bresx * bresy;
+ G_debug(2, " barea,sarea,rwalk,sisum: %f %f %f %f", barea, sarea,
+ rwalk, sisum);
+ /* write hh.walkers0 */
+#ifdef PARALLEL
+ G_message(_("PREPARE"));
+#endif
+
+
+ G_message(_("my %d,mx %d"),my,mx);
+/* for (cnt1 = 0; cnt1 < my; cnt1++) {
+ for (cnt2 = 0; cnt2 < mx; cnt2++) {
+
+
+ //G_message(_("gen"));
+ gen = rwalk * si[cnt1][cnt2] / sisum;
+ // G_message(_("mgen"));
+ mgen = (int)gen;
+ // G_message(_("if"));
+ if (cnt2 == 0)
+ {
+ if (cnt1 == 0)
+ {
+ // G_message(_("ifif"));
+ lw_start[cnt1][cnt2] = (zz[cnt1][cnt2] != UNDEF ? mgen : 0);
+ }
+ else
+ {
+ // G_message(_("ifelse"));
+ lw_start[cnt1][cnt2] = lw_start[cnt1 - 1][mx - 1] + (zz[cnt1][cnt2] != UNDEF ? mgen : 0);
+ }
+ }
+ else
+ {
+ // G_message(_("ifelse"));
+ lw_start[cnt1][cnt2] = lw_start[cnt1][cnt2 - 1] + (zz[cnt1][cnt2] != UNDEF ? mgen : 0);
+ }
+
+ // G_message(_("next"));
+ }
+ }*/
+ /* write hh.walkers0 */
+#ifdef PARALLEL
+ G_message(_("PARAWORK1"));
+#endif
+
+//---------------------------------------------------------------------------------------------------------------
+#ifdef PARALLEL
+ printTimeDiff("P0");
+#endif
+ /* write hh.walkers0 */
+ //recalculate: lw
+ // #pragma omp parallel firstprivate(gen,mgen,wei,x,y,gen2,mgen2,wei2,mgen3,wei3,nmult,weifac,iw,walkwe)
+
+//{
+ // #pragma omp for schedule(dynamic) reduction(max:lw)
+ for (k = 0; k < my; k++) {
+ for (l = 0; l < mx; l++) { /* run thru the whole area */
+ if (zz[k][l] != UNDEF) {
+
+ x = xp0 + stepx * (double)(l);
+ y = yp0 + stepy * (double)(k);
+
+ gen = rwalk * si[k][l] / sisum;
+ mgen = (int)gen;
+ wei = gen / (double)(mgen + 1);
+
+ /*if (si[k][l] != 0.) { */
+ /* this stuff later for multiscale */
+
+ gen2 =
+ (double)maxwab *si[k][l] / (si0 *
+ (double)(mx2o * my2o));
+ gen2 = gen2 * (barea / sarea);
+ mgen2 = (int)gen2;
+ wei2 = gen2 / (double)(mgen2 + 1);
+ mgen3 =
+ (int)((double)mgen2 * wei2 / ((double)mgen * wei));
+ nmult = mgen3 + 1;
+ wei3 = gen2 / (double)((mgen + 1) * (mgen2 + 1));
+ weifac = wei3 / wei;
+ /* } else {
+ nmult = 1;
+ weifac = 1.;
+ fprintf(stderr, "\n zero rainfall excess in cell");
+ } */
+
+ /*G_debug(2, " gen,gen2,wei,wei2,mgen3,nmult: %f %f %f %f %d %d",gen,gen2,wei,wei2,mgen3,nmult);
+ */
+ // lw = lw_start[k][l];
+ //if (lw != lw_start[k][l])
+ //
+ // G_message(_("lw %d lws %d k %d l %d"),lw,lw_start[k][l],k,l);
+ //}
+ for (iw = 1; iw <= mgen + 1; iw++) { /* assign walkers */
+
+ if (lw >= MAXW) /* max valid value is MAXW - 1, not MAXW */
+ G_fatal_error(_("nwalk (%d) > maxw (%d)!"), lw, MAXW);
+
+ w[lw][0] = x + stepx * (simwe_rand() - 0.5);
+ w[lw][1] = y + stepy * (simwe_rand() - 0.5);
+ w[lw][2] = wei;
+
+ walkwe += w[lw][2];
+ vavg[lw][0] = v1[k][l];
+ vavg[lw][1] = v2[k][l];
+ if (w[lw][0] >= xmin && w[lw][1] >= ymin &&
+ w[lw][0] <= xmax && w[lw][1] <= ymax) {
+ iflag[lw] = 0;
+ }
+ else {
+ iflag[lw] = 1;
+ }
+ lw++;
+ }
+ } /*DEFined area */
+ }
+ }
+//}
+ nwalk = lw;
+ G_debug(2, " nwalk, maxw %d %d", nwalk, MAXW);
+ G_debug(2, " walkwe (walk weight),frac %f %f", walkwe, frac);
+#ifdef PARALLEL
+ printTimeDiff("P1");
+#endif
+#ifdef PARALLEL
+ G_message(_("PARAWORK2"));
+#endif
+ stxm = stepx * (double)(mx + 1) - xmin;
+ stym = stepy * (double)(my + 1) - ymin;
+ nwalka = 0;
+ deldif = sqrt(deltap) * frac; /* diffuse factor */
+
+
+ factor = deltap * sisum / (rwalk * (double)nblock);
+
+ G_debug(2, " deldif,factor %f %e", deldif, factor);
+
+ /* ********************************************************** */
+ /* main loop over the projection time */
+ /* *********************************************************** */
+
+
+ G_debug(2, "main loop over the projection time... ");
+
+
+ //double *gaux_arr = (double *)G_calloc((nwalk), sizeof(double *));;
+ //double *gauy_arr = (double *)G_calloc((nwalk), sizeof(double *));;
+ // int stop = -1;
+G_message("miter %d",miter);
+ for (i = 1; i <= miter; i++) { /* iteration loop depending on simulation time and deltap */
+ /* #pragma omp critical STOP
+{
+ if (stop != -1 && i >= stop)
+ {
+ continue;
+ }
+}*/
+
+ G_percent(i, miter, 1);
+ iter1 = i / iterout;
+ iter1 *= iterout;
+ if (iter1 == i) {
+ nfiterw = i / iterout + 10;
+ nfiterh = i / iterout + 40;
+ G_debug(2, "iblock=%d i=%d miter=%d nwalk=%d nwalka=%d",
+ iblock, i, miter, nwalk, nwalka);
+ }
+
+ if (nwalka == 0 && i > 1)
+ {
+ /* #pragma omp critical STOP
+{
+ if (stop != -1 && i >= stop)
+ {
+ continue;
+ }
+}*/
+ goto L_800;
+ }
+ /* ************************************************************ */
+ /* .... propagate one step */
+ /* ************************************************************ */
+
+ addac = factor;
+ conn = (double)nblock / (double)iblock;
+ if (i == 1) {
+ addac = factor * .5;
+ }
+ nwalka = 0;
+ nstack = 0;
+
+ /*for (lw = 0; lw < nwalk; lw++) {
+ if (w[lw][2] > EPS) {
+ ++nwalka;
+ gaux_arr[lw] = gasdev();
+ gauy_arr[lw] = gasdev();
+ }
+ }*/
+ //int chunk_size = (int)nwalk/(4 * omp_get_thread_num());
+
+
+
+#pragma omp parallel firstprivate(l,lw,k,decr,d1,hhc,velx,vely,eff,gaux,gauy)//nwalka
+{
+ int steps = (int)((((double)nwalk) / ((double) omp_get_num_threads())) + 0.5);
+ int tid = omp_get_thread_num();
+ int min_loop = tid * steps;
+ int max_loop = ((tid + 1) * steps) > nwalk ? nwalk : (tid + 1) * steps;
+/*#pragma omp critical
+ {
+ G_message("tid: %d min_loop: %d max_loop: %d nwalk: %d, steps: %d", tid, min_loop, max_loop, nwalk, steps);
+ }*/
+ //#pragma omp for schedule(static,1)
+ for (lw = min_loop; lw < max_loop; lw++) {
+ //G_message(_("nthreads %d"),omp_get_num_threads());
+ //G_message("\n tid: %d", omp_get_thread_num());
+ if (w[lw][2] > EPS) { /* check the walker weight */
+ //#pragma omp atomic
+ ++nwalka;
+ l = (int)((w[lw][0] + stxm) / stepx) - mx - 1;
+ k = (int)((w[lw][1] + stym) / stepy) - my - 1;
+//G_message("k_%d l_%d", k,l);
+ if (l > mx - 1 || k > my - 1 || k < 0 || l < 0) {
+
+ G_debug(2, " k,l=%d,%d", k, l);
+ printf(" lw,w=%d %f %f", lw, w[lw][1], w[lw][2]);
+ G_debug(2, " stxym=%f %f", stxm, stym);
+ printf(" step=%f %f", stepx, stepy);
+ G_debug(2, " m=%d %d", my, mx);
+ printf(" nwalka,nwalk=%d %d", nwalka, nwalk);
+ G_debug(2, " ");
+ }
+
+ if (zz[k][l] != UNDEF) {
+ if (infil != NULL) { /* infiltration part */
+ if (inf[k][l] - si[k][l] > 0.) {
+
+ decr = pow(addac * w[lw][2], 3. / 5.); /* decreasing factor in m */
+// #pragma omp critical
+//{
+ if (inf[k][l] > decr) {
+ inf[k][l] -= decr; /* decrease infilt. in cell and eliminate the walker */
+ w[lw][2] = 0.;
+ }
+ else {
+ w[lw][2] -= pow(inf[k][l], 5. / 3.) / addac; /* use just proportional part of the walker weight */
+ inf[k][l] = 0.;
+
+ }
+//}
+ }
+ }
+
+//#pragma omp critical
+//{
+ gama[k][l] += (addac * w[lw][2]); /* add walker weigh to water depth or conc. */
+ d1 = gama[k][l] * conn;
+//
+//}
+ gasdev_for_paralel(&gaux,&gauy);
+ hhc = pow(d1, 3. / 5.);
+
+ if (hhc > hhmax && wdepth == NULL) { /* increased diffusion if w.depth > hhmax */
+ dif[k][l] = (halpha + 1) * deldif;
+ velx = vavg[lw][0];
+ vely = vavg[lw][1];
+ }
+ else {
+ dif[k][l] = deldif;
+ velx = v1[k][l];
+ vely = v2[k][l];
+ }
+
+
+ if (traps != NULL && trap[k][l] != 0.) { /* traps */
+
+ eff = simwe_rand(); /* random generator */
+
+ if (eff <= trap[k][l]) {
+ velx = -0.1 * v1[k][l]; /* move it slightly back */
+ vely = -0.1 * v2[k][l];
+ }
+ }
+
+
+ w[lw][0] += (velx + dif[k][l] * gaux); /* move the walker */
+ w[lw][1] += (vely + dif[k][l] * gauy);
+
+ if (hhc > hhmax && wdepth == NULL) {
+ vavg[lw][0] = hbeta * (vavg[lw][0] + v1[k][l]);
+ vavg[lw][1] = hbeta * (vavg[lw][1] + v2[k][l]);
+ }
+
+ if (w[lw][0] <= xmin || w[lw][1] <= ymin || w[lw][0]
+ >= xmax || w[lw][1] >= ymax) {
+ w[lw][2] = 1e-10; /* eliminate walker if it is out of area */
+ }
+ else {
+ if (wdepth != NULL) {
+ l = (int)((w[lw][0] + stxm) / stepx) - mx - 1;
+ k = (int)((w[lw][1] + stym) / stepy) - my - 1;
+ w[lw][2] *= sigma[k][l];
+ }
+
+ } /* else */
+ } /*DEFined area */
+ else {
+ w[lw][2] = 1e-10; /* eliminate walker if it is out of area */
+ }
+ }
+ } /* lw loop */
+}
+ /* Changes made by Soeren 8. Mar 2011 to replace the site walker output implementation */
+ /* Save all walkers located within the computational region and with valid
+ z coordinates */
+ if ((i == miter || i == iter1)) {
+ nstack = 0;
+
+ for (lw = 0; lw < nwalk; lw++) {
+ /* Compute the elevation raster map index */
+ l = (int)((w[lw][0] + stxm) / stepx) - mx - 1;
+ k = (int)((w[lw][1] + stym) / stepy) - my - 1;
+
+ /* Check for correct elevation raster map index */
+ if(l < 0 || l >= mx || k < 0 || k >= my)
+ continue;
+
+ if (w[lw][2] > EPS && zz[k][l] != UNDEF) {
+
+ /* Save the 3d position of the walker */
+ stack[nstack][0] = mixx / conv + w[lw][0] / conv;
+ stack[nstack][1] = miyy / conv + w[lw][1] / conv;
+ stack[nstack][2] = zz[k][l];
+
+ nstack++;
+ }
+ } /* lw loop */
+ }
+
+ if (i == iter1 && ts == 1) {
+ /* call output for iteration output */
+ if (erdep != NULL)
+ erod(gama); /* divergence of gama field */
+
+ conn = (double)nblock / (double)iblock;
+ itime = (int)(i * deltap * timec);
+ ii = output_data(itime, conn);
+ if (ii != 1)
+ G_fatal_error(_("Unable to write raster maps"));
+ }
+
+ /* Write the water depth each time step at an observation point */
+ if(points.is_open)
+ {
+ double value = 0.0;
+ int p;
+ fprintf(points.output, "%.6d ", i);
+ /* Write for each point */
+ for(p = 0; p < points.npoints; p++)
+ {
+ l = (int)((points.x[p] - mixx + stxm) / stepx) - mx - 1;
+ k = (int)((points.y[p] - miyy + stym) / stepy) - my - 1;
+
+ if (zz[k][l] != UNDEF) {
+
+ if (wdepth == NULL)
+ value = step * gama[k][l] * cchez[k][l];
+ else
+ value = gama[k][l] * slope[k][l];
+
+ fprintf(points.output, "%2.4f ", value);
+ } else {
+ /* Point is invalid, so a negative value is written */
+ fprintf(points.output, "%2.4f ", -1.0);
+ }
+ }
+ fprintf(points.output, "\n");
+ }
+ } /* miter */
+ L_800:
+ /* Soeren 8. Mar 2011: Why is this commented out?*/
+ /* if (iwrib != nblock) {
+ icount = icoub / iwrib;
+
+ if (icoub == (icount * iwrib)) {
+ ++icfl;
+ nflw = icfl + 50;
+ conn = (double) nblock / (double) iblock;
+
+ }
+ } */
+
+
+#ifdef PARALLEL
+ G_message(_("PARAWORK3"));
+#endif
+#ifdef PARALLEL
+ printTimeDiff("P2");
+#endif
+ if (err != NULL) {
+ for (k = 0; k < my; k++) {
+ for (l = 0; l < mx; l++) {
+ if (zz[k][l] != UNDEF) {
+ d1 = gama[k][l] * (double)conn;
+ gammas[k][l] += pow(d1, 3. / 5.);
+ } /* DEFined area */
+ }
+ }
+ }
+ if (erdep != NULL)
+ erod(gama);
+ }
+ /* ........ end of iblock loop */
+
+#ifdef PARALLEL
+ printTimeDiff("L1");
+#endif
+ /* Write final maps here because we know the last time stamp here */
+ if (ts == 0) {
+ conn = (double)nblock / (double)iblock;
+ itime = (int)(i * deltap * timec);
+ ii = output_data(itime, conn);
+ if (ii != 1)
+ G_fatal_error(_("Cannot write raster maps"));
+ }
+ /* Close the observation logfile */
+ if(points.is_open)
+ fclose(points.output);
+#ifdef PARALLEL
+ printTimeDiff("L1");
+#endif
+ points.is_open = 0;
+
+}
Added: grass-addons/grass7/raster/r.sim.water.mp/input.c
===================================================================
--- grass-addons/grass7/raster/r.sim.water.mp/input.c (rev 0)
+++ grass-addons/grass7/raster/r.sim.water.mp/input.c 2017-04-28 04:43:41 UTC (rev 70973)
@@ -0,0 +1,747 @@
+/* input.c (simlib), 20.nov.2002, JH */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include <grass/linkm.h>
+#include <grass/gmath.h>
+#include "simlib.h"
+#include "waterglobs.h"
+
+
+/* Local prototypes for raster map reading and array allocation */
+static float ** read_float_raster_map(int rows, int cols, char *name, float unitconv);
+static double ** read_double_raster_map(int rows, int cols, char *name, double unitconv);
+static float ** create_float_matrix(int rows, int cols, float fill_value);
+static double ** create_double_matrix(int rows, int cols, double fill_value);
+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 */
+/* *************************************************************** */
+
+/*!
+ * \brief allocate memory, read input rasters, assign UNDEF to NODATA
+ *
+ * \return int
+ */
+
+/* ************************************************************************* */
+/* Read all input maps and input values into memory ************************ */
+int input_data(void)
+{
+ int rows = my, cols = mx; /* my and mx are global variables */
+ double unitconv = 0.000000278; /* mm/hr to m/s */
+ int if_rain = 0;
+
+ G_debug(1, "Running MAR 2011 version, started modifications on 20080211");
+ G_debug(1, "Reading input data");
+
+ /* Elevation and gradients are mandatory */
+ zz = read_float_raster_map(rows, cols, elevin, 1.0);
+ v1 = read_double_raster_map(rows, cols, dxin, 1.0);
+ v2 = read_double_raster_map(rows, cols, dyin, 1.0);
+
+ /* Update elevation map */
+ copy_matrix_undef_double_to_float_values(rows, cols, v1, zz);
+ copy_matrix_undef_double_to_float_values(rows, cols, v2, zz);
+
+ /* Manning surface roughnes: read map or use a single value */
+ if(manin != NULL) {
+ cchez = read_float_raster_map(rows, cols, manin, 1.0);
+ } else if(manin_val >= 0.0) { /* If no value set its set to -999.99 */
+ cchez = create_float_matrix(rows, cols, manin_val);
+ }else{
+ G_fatal_error(_("Raster map <%s> not found, and manin_val undefined, choose one to be allowed to process"), manin);
+ }
+
+ /* Rain: read rain map or use a single value for all cells */
+ if (rain != NULL) {
+ si = read_double_raster_map(rows, cols, rain, unitconv);
+ if_rain = 1;
+ } else if(rain_val >= 0.0) { /* If no value set its set to -999.99 */
+ si = create_double_matrix(rows, cols, rain_val * unitconv);
+ if_rain = 1;
+ } else{
+ si = create_double_matrix(rows, cols, (double)UNDEF);
+ if_rain = 0;
+ }
+
+ /* Update elevation map */
+ copy_matrix_undef_double_to_float_values(rows, cols, si, zz);
+
+ /* Load infiltration and traps if rain is present */
+ if(if_rain == 1) {
+ /* Infiltration: read map or use a single value */
+ if (infil != NULL) {
+ inf = read_double_raster_map(rows, cols, infil, unitconv);
+ } else if(infil_val >= 0.0) { /* If no value set its set to -999.99 */
+ inf = create_double_matrix(rows, cols, infil_val * unitconv);
+ } else{
+ inf = create_double_matrix(rows, cols, (double)UNDEF);
+ }
+
+ /* Traps */
+ if (traps != NULL)
+ trap = read_float_raster_map(rows, cols, traps, 1.0);
+ else
+ trap = create_float_matrix(rows, cols, (double)UNDEF);
+ }
+
+ if (detin != NULL) {
+ dc = read_float_raster_map(rows, cols, detin, 1.0);
+ copy_matrix_undef_float_values(rows, cols, dc, zz);
+ }
+
+ if (tranin != NULL) {
+ ct = read_float_raster_map(rows, cols, tranin, 1.0);
+ copy_matrix_undef_float_values(rows, cols, ct, zz);
+ }
+
+ if (tauin != NULL) {
+ tau = read_float_raster_map(rows, cols, tauin, 1.0);
+ copy_matrix_undef_float_values(rows, cols, tau, zz);
+ }
+
+ if (wdepth != NULL) {
+ gama = read_double_raster_map(rows, cols, wdepth, 1.0);
+ copy_matrix_undef_double_to_float_values(rows, cols, gama, zz);
+ }
+
+ /* Array for gradient checking */
+ slope = create_double_matrix(rows, cols, 0.0);
+
+ /* Create the observation points and open the logfile */
+ create_observation_points();
+
+ return 1;
+}
+
+/* ************************************************************************* */
+
+/* data preparations, sigma, shear, etc. */
+int grad_check(void)
+{
+ int k, l;
+ double zx, zy, zd2, zd4, sinsl;
+ double cc, cmul2;
+ double sheer;
+ double vsum = 0.;
+ double vmax = 0.;
+ double chsum = 0.;
+ double zmin = 1.e12;
+ double zmax = -1.e12;
+ double zd2min = 1.e12;
+ double zd2max = -1.e12;
+ double smin = 1.e12;
+ double smax = -1.e12;
+ double infmin = 1.e12;
+ double infmax = -1.e12;
+ double sigmax = -1.e12;
+ double cchezmax = -1.e12;
+ double rhow = 1000.;
+ double gacc = 9.81;
+ double hh = 1.;
+ double deltaw = 1.e12;
+
+ sisum = 0.;
+ infsum = 0.;
+ cmul2 = rhow * gacc;
+
+ for (k = 0; k < my; k++) {
+ for (l = 0; l < mx; l++) {
+ if (zz[k][l] != UNDEF) {
+ zx = v1[k][l];
+ zy = v2[k][l];
+ zd2 = zx * zx + zy * zy;
+ sinsl = sqrt(zd2) / sqrt(zd2 + 1); /* sin(terrain slope) */
+ /* Computing MIN */
+ zd2 = sqrt(zd2);
+ zd2min = amin1(zd2min, zd2);
+ /* Computing MAX */
+ zd2max = amax1(zd2max, zd2);
+ zd4 = sqrt(zd2); /* ^.25 */
+ if (cchez[k][l] != 0.) {
+ cchez[k][l] = 1. / cchez[k][l]; /* 1/n */
+ }
+ else {
+ G_fatal_error(_("Zero value in Mannings n"));
+ }
+ if (zd2 == 0.) {
+ v1[k][l] = 0.;
+ v2[k][l] = 0.;
+ slope[k][l] = 0.;
+ }
+ else {
+ if (wdepth)
+ hh = pow(gama[k][l], 2. / 3.);
+ /* hh = 1 if there is no water depth input */
+ v1[k][l] = (double)hh *cchez[k][l] * zx / zd4;
+ v2[k][l] = (double)hh *cchez[k][l] * zy / zd4;
+
+ slope[k][l] =
+ sqrt(v1[k][l] * v1[k][l] + v2[k][l] * v2[k][l]);
+ }
+ if (wdepth) {
+ sheer = (double)(cmul2 * gama[k][l] * sinsl); /* shear stress */
+ /* if critical shear stress >= shear then all zero */
+ if ((sheer <= tau[k][l]) || (ct[k][l] == 0.)) {
+ si[k][l] = 0.;
+ sigma[k][l] = 0.;
+ }
+ else {
+ si[k][l] = (double)(dc[k][l] * (sheer - tau[k][l]));
+ sigma[k][l] = (double)(dc[k][l] / ct[k][l]) * (sheer - tau[k][l]) / (pow(sheer, 1.5)); /* rill erosion=1.5, sheet = 1.1 */
+ }
+ }
+ sisum += si[k][l];
+ smin = amin1(smin, si[k][l]);
+ smax = amax1(smax, si[k][l]);
+ if (inf) {
+ infsum += inf[k][l];
+ infmin = amin1(infmin, inf[k][l]);
+ infmax = amax1(infmax, inf[k][l]);
+ }
+ vmax = amax1(vmax, slope[k][l]);
+ vsum += slope[k][l];
+ chsum += cchez[k][l];
+ zmin = amin1(zmin, (double)zz[k][l]);
+ zmax = amax1(zmax, (double)zz[k][l]); /* not clear were needed */
+ if (wdepth)
+ sigmax = amax1(sigmax, sigma[k][l]);
+ cchezmax = amax1(cchezmax, cchez[k][l]);
+ /* saved sqrt(sinsl)*cchez to cchez array for output */
+ cchez[k][l] *= sqrt(sinsl);
+ } /* DEFined area */
+ }
+ }
+ if (inf != NULL && smax < infmax)
+ G_warning(_("Infiltration exceeds the rainfall rate everywhere! No overland flow."));
+
+ cc = (double)mx *my;
+
+ si0 = sisum / cc;
+ vmean = vsum / cc;
+ chmean = chsum / cc;
+
+ if (inf)
+ infmean = infsum / cc;
+
+ if (wdepth)
+ deltaw = 0.8 / (sigmax * vmax); /*time step for sediment */
+ deltap = 0.25 * sqrt(stepx * stepy) / vmean; /*time step for water */
+
+ if (deltaw > deltap)
+ timec = 4.;
+ else
+ timec = 1.25;
+
+ miter = (int)(timesec / (deltap * timec)); /* number of iterations = number of cells to pass */
+ iterout = (int)(iterout / (deltap * timec)); /* number of cells to pass for time series output */
+
+ fprintf(stderr, "\n");
+ G_message(_("Min elevation \t= %.2f m\nMax elevation \t= %.2f m\n"), zmin,
+ zmax);
+ G_message(_("Mean Source Rate (rainf. excess or sediment) \t= %f m/s or kg/m2s \n"),
+ si0);
+ G_message(_("Mean flow velocity \t= %f m/s\n"), vmean);
+ G_message(_("Mean Mannings \t= %f\n"), 1.0 / chmean);
+
+ deltap = amin1(deltap, deltaw);
+
+ G_message(n_("Number of iterations \t= %d cell\n",
+ "Number of iterations \t= %d cells\n", miter), miter);
+ G_message(_("Time step \t= %.2f s\n"), deltap);
+ if (wdepth) {
+ G_message(_("Sigmax \t= %f\nMax velocity \t= %f m/s\n"), sigmax,
+ vmax);
+ G_message(_("Time step used \t= %.2f s\n"), deltaw);
+ }
+ /* if (wdepth) deltap = 0.1;
+ * deltap for sediment is ar. average deltap and deltaw */
+ /* if (wdepth) deltap = (deltaw+deltap)/2.;
+ * deltap for sediment is ar. average deltap and deltaw */
+
+
+ /*! For each cell (k,l) compute the length s=(v1,v2) of the path
+ * that the particle will travel per one time step
+ * \f$ s(k,l)=v(k,l)*dt \f$, [m]=[m/s]*[s]
+ * give warning if there is a cell that will lead to path longer than 2 cells
+ *
+ * if running erosion, compute sediment transport capacity for each cell si(k,l)
+ * \f$
+ * T({\bf r})=K_t({\bf r}) \bigl[\tau({\bf r})\bigr]^p
+ * =K_t({\bf r}) \bigl[\rho_w\, g h({\bf r}) \sin \beta ({\bf r}) \bigr]^p
+ * \f$
+ * [kg/ms]=...
+ */
+ for (k = 0; k < my; k++) {
+ for (l = 0; l < mx; l++) {
+ if (zz[k][l] != UNDEF) {
+ v1[k][l] *= deltap;
+ v2[k][l] *= deltap;
+ /*if(v1[k][l]*v1[k][l]+v2[k][l]*v2[k][l] > cellsize, warning, napocitaj
+ *ak viac ako 10%a*/
+ /* THIS IS CORRECT SOLUTION currently commented out */
+ if (inf)
+ inf[k][l] *= timesec;
+ if (wdepth)
+ gama[k][l] = 0.;
+ if (et) {
+ if (sigma[k][l] == 0. || slope[k][l] == 0.)
+ si[k][l] = 0.;
+ else
+ /* temp for transp. cap. erod */
+ si[k][l] = si[k][l] / (slope[k][l] * sigma[k][l]);
+ }
+ } /* DEFined area */
+ }
+ }
+
+ /*! compute transport capacity limted erosion/deposition et
+ * as a divergence of sediment transport capacity
+ * \f$
+ D_T({\bf r})= \nabla\cdot {\bf T}({\bf r})
+ * \f$
+ */
+ if (et) {
+ erod(si); /* compute divergence of t.capc */
+ if (output_et() != 1)
+ G_fatal_error(_("Unable to write et file"));
+ }
+
+ /*! compute the inversion operator and store it in sigma - note that after this
+ * sigma does not store the first order reaction coefficient but the operator
+ * WRITE the equation here
+ */
+ if (wdepth) {
+ for (k = 0; k < my; k++) {
+ for (l = 0; l < mx; l++) {
+ if (zz[k][l] != UNDEF) {
+ /* get back from temp */
+ if (et)
+ si[k][l] = si[k][l] * slope[k][l] * sigma[k][l];
+ if (sigma[k][l] != 0.)
+ /* rate of weight loss - w=w*sigma ,
+ * vaha prechadzky po n-krokoch je sigma^n */
+
+ /*!!!!! not clear what's here :-\ !!!!!*/
+
+ sigma[k][l] =
+ exp(-sigma[k][l] * deltap * slope[k][l]);
+ /* if(sigma[k][l]<0.5) warning, napocitaj,
+ * ak vacsie ako 50% skonci, zmensi deltap)*/
+ }
+ } /*DEFined area */
+ }
+ }
+ return 1;
+}
+
+/* ************************************************************************* */
+
+void copy_matrix_undef_double_to_float_values(int rows, int cols, double **source, float **target)
+{
+ int col = 0, row = 0;
+
+ for(row = 0; row < rows; row++) {
+ for(col = 0; col < cols; col++) {
+ if(source[row][col] == UNDEF)
+ target[row][col] = UNDEF;
+ }
+ }
+}
+
+/* ************************************************************************* */
+
+void copy_matrix_undef_float_values(int rows, int cols, float **source, float **target)
+{
+ int col = 0, row = 0;
+
+ for(row = 0; row < rows; row++) {
+ for(col = 0; col < cols; col++) {
+ if(source[row][col] == UNDEF)
+ target[row][col] = UNDEF;
+ }
+ }
+}
+
+/* ************************************************************************* */
+
+float ** create_float_matrix(int rows, int cols, float fill_value)
+{
+ int col = 0, row = 0;
+ float **matrix = NULL;
+
+ G_verbose_message("Creating float matrix with value %g", fill_value);
+
+ /* Allocate the float marix */
+ matrix = G_alloc_fmatrix(rows, cols);
+
+ for(row = 0; row < rows; row++) {
+ for(col = 0; col < cols; col++) {
+ matrix[row][col] = fill_value;
+ }
+ }
+
+ return matrix;
+}
+
+/* ************************************************************************* */
+
+double ** create_double_matrix(int rows, int cols, double fill_value)
+{
+ int col = 0, row = 0;
+ double **matrix = NULL;
+
+ G_verbose_message("Creating double matrix with value %g", fill_value);
+
+ /* Allocate the float marix */
+ matrix = G_alloc_matrix(rows, cols);
+
+ for(row = 0; row < rows; row++) {
+ for(col = 0; col < cols; col++) {
+ matrix[row][col] = fill_value;
+ }
+ }
+
+ return matrix;
+}
+
+/* ************************************************************************* */
+
+float ** read_float_raster_map(int rows, int cols, char *name, float unitconv)
+{
+ FCELL *row_buff = NULL;
+ int fd;
+ int col = 0, row = 0, row_rev = 0;
+ float **matrix = NULL;
+
+ G_verbose_message("Reading float map %s into memory", name);
+
+ /* Open raster map */
+ fd = Rast_open_old(name, "");
+
+ /* Allocate the row buffer */
+ row_buff = Rast_allocate_f_buf();
+
+ /* Allocate the float marix */
+ matrix = G_alloc_fmatrix(rows, cols);
+
+ for(row = 0; row < rows; row++) {
+ Rast_get_f_row(fd, row_buff, row);
+
+ for(col = 0; col < cols; col++) {
+ /* we fill the arrays from south to north */
+ row_rev = rows - row - 1;
+ /* Check for null values */
+ if (!Rast_is_f_null_value(row_buff + col))
+ matrix[row_rev][col] = (float)(unitconv * row_buff[col]);
+ else
+ matrix[row_rev][col] = UNDEF;
+ }
+ }
+
+ /* Free the row buffer */
+ if(row_buff)
+ G_free(row_buff);
+
+ Rast_close(fd);
+
+ return matrix;
+}
+
+/* ************************************************************************* */
+
+double ** read_double_raster_map(int rows, int cols, char *name, double unitconv)
+{
+ DCELL *row_buff = NULL;
+ int fd;
+ int col = 0, row = 0, row_rev;
+ double **matrix = NULL;
+
+ G_verbose_message("Reading double map %s into memory", name);
+
+ /* Open raster map */
+ fd = Rast_open_old(name, "");
+
+ /* Allocate the row buffer */
+ row_buff = Rast_allocate_d_buf();
+
+ /* Allocate the double marix */
+ matrix = G_alloc_matrix(rows, cols);
+
+ for(row = 0; row < rows; row++) {
+ Rast_get_d_row(fd, row_buff, row);
+
+ for(col = 0; col < cols; col++) {
+ /* we fill the arrays from south to north */
+ row_rev = rows - row - 1;
+ /* Check for null values */
+ if (!Rast_is_d_null_value(row_buff + col))
+ matrix[row_rev][col] = (double)(unitconv * row_buff[col]);
+ else
+ matrix[row_rev][col] = UNDEF;
+ }
+ }
+
+ /* Free the row buffer */
+ if(row_buff)
+ G_free(row_buff);
+
+ Rast_close(fd);
+
+ return matrix;
+}
Added: grass-addons/grass7/raster/r.sim.water.mp/main.c
===================================================================
--- grass-addons/grass7/raster/r.sim.water.mp/main.c (rev 0)
+++ grass-addons/grass7/raster/r.sim.water.mp/main.c 2017-04-28 04:43:41 UTC (rev 70973)
@@ -0,0 +1,538 @@
+
+/****************************************************************************
+ *
+ * MODULE: r.sim.water: main program for hydrologic and sediment transport
+ * simulation (SIMWE)
+ *
+ * AUTHOR(S): L. Mitas, H. Mitasova, J. Hofierka
+ * PURPOSE: Hydrologic and sediment transport simulation (SIMWE)
+ *
+ * COPYRIGHT: (C) 2002, 2010 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.
+ *
+ *****************************************************************************/
+
+/*-
+ * r.sim.water.mp: main program for hydrologic and sediment transport
+ * simulation (SIMWE)
+ *
+ * Original program (2002) and various modifications:
+ * Lubos Mitas, Helena Mitasova
+ *
+ * GRASS5.0 version of the program:
+ * J. Hofierka
+ *
+ * Copyright (C) 2002 L. Mitas, H. Mitasova, J. Hofierka
+ *
+ *This program is free software; you can redistribute it and/or
+ *modify it under the terms of the GNU General Public License
+ *as published by the Free Software Foundation; either version 2
+ *of the License, or (at your option) any later version.
+ *
+ *This program is distributed in the hope that it will be useful,
+ *but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ *GNU General Public License for more details.
+ *
+ *You should have received a copy of the GNU General Public License
+ *along with this program; if not, write to the Free Software
+ *Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+ *
+ *
+ * Notes on modifications:
+ * v. 1.0 May 2002
+ * modified by Y. Chemin in February 2008 (reporting, optional inputs)
+ * sites-related input/output commented out Nov. 2008
+ * December 2016, OpenMP support Stanislav Zubal, Michal Lacko
+ */
+
+/********************************/
+/* DEFINE GLOB VAR */
+
+/********************************/
+/* #define NWALK "1000000" */
+#define DIFFC "0.8"
+#define HMAX "0.3"
+#define HALPHA "4.0"
+#define HBETA "0.5"
+#define NITER "10"
+#define ITEROUT "2"
+#define DENSITY "200"
+#define RAINVAL "50"
+#define MANINVAL "0.1"
+#define INFILVAL "0.0"
+
+/********************************/
+/* INCLUDES */
+
+/********************************/
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+#include <grass/config.h>
+#ifdef HAVE_UNISTD_H
+#include <unistd.h>
+#endif
+#include <grass/gis.h>
+#include <grass/vector.h>
+#include <grass/linkm.h>
+#include <grass/bitmap.h>
+#include <grass/glocale.h>
+#include <grass/gmath.h>
+
+/********************************/
+/* Specific stuff */
+
+/********************************/
+
+#include "simlib.h"
+
+
+/****************************************/
+/* MAIN */
+
+/****************************************/
+
+#ifdef PARALLEL
+
+void printTimeDiff(const char* message)
+{
+ double seconds;
+ time_t timer2;
+ time(&timer2); //or timer2 = time(NULL);
+ seconds = difftime(timer2,timer);
+ timer = timer2;
+ G_message(_("TIME_DIFF: <%.0lf>\tmessage <%s>\n"),
+ seconds,
+ message);
+}
+#endif
+
+int main(int argc, char *argv[])
+{
+ int ii;
+ int ret_val;
+ 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]);
+
+
+#ifdef PARALLEL
+ timer = time(NULL);
+ printTimeDiff("Start");
+#endif
+
+ module = G_define_module();
+ G_add_keyword(_("raster"));
+ G_add_keyword(_("hydrology"));
+ G_add_keyword(_("soil"));
+ G_add_keyword(_("flow"));
+ G_add_keyword(_("overland flow"));
+ G_add_keyword(_("model"));
+ module->description =
+ _("Overland flow hydrologic simulation using "
+ "path sampling method (SIMWE).");
+
+ parm.elevin = G_define_standard_option(G_OPT_R_ELEV);
+
+ parm.dxin = G_define_standard_option(G_OPT_R_INPUT);
+ parm.dxin->key = "dx";
+ parm.dxin->description = _("Name of x-derivatives raster map [m/m]");
+
+ parm.dyin = G_define_standard_option(G_OPT_R_INPUT);
+ parm.dyin->key = "dy";
+ parm.dyin->description = _("Name of y-derivatives raster map [m/m]");
+
+ parm.rain = G_define_standard_option(G_OPT_R_INPUT);
+ parm.rain->key = "rain";
+ parm.rain->required = NO;
+ parm.rain->description =
+ _("Name of rainfall excess rate (rain-infilt) raster map [mm/hr]");
+ parm.rain->guisection = _("Input");
+
+ parm.rainval = G_define_option();
+ parm.rainval->key = "rain_value";
+ parm.rainval->type = TYPE_DOUBLE;
+ parm.rainval->answer = RAINVAL;
+ parm.rainval->required = NO;
+ parm.rainval->description =
+ _("Rainfall excess rate unique value [mm/hr]");
+ parm.rainval->guisection = _("Input");
+
+ parm.infil = G_define_standard_option(G_OPT_R_INPUT);
+ parm.infil->key = "infil";
+ parm.infil->required = NO;
+ parm.infil->description =
+ _("Name of runoff infiltration rate raster map [mm/hr]");
+ parm.infil->guisection = _("Input");
+
+ parm.infilval = G_define_option();
+ parm.infilval->key = "infil_value";
+ parm.infilval->type = TYPE_DOUBLE;
+ parm.infilval->answer = INFILVAL;
+ parm.infilval->required = NO;
+ parm.infilval->description =
+ _("Runoff infiltration rate unique value [mm/hr]");
+ parm.infilval->guisection = _("Input");
+
+ parm.manin = G_define_standard_option(G_OPT_R_INPUT);
+ parm.manin->key = "man";
+ parm.manin->required = NO;
+ parm.manin->description = _("Name of Manning's n raster map");
+ parm.manin->guisection = _("Input");
+
+ parm.maninval = G_define_option();
+ parm.maninval->key = "man_value";
+ parm.maninval->type = TYPE_DOUBLE;
+ parm.maninval->answer = MANINVAL;
+ parm.maninval->required = NO;
+ parm.maninval->description = _("Manning's n unique value");
+ parm.maninval->guisection = _("Input");
+
+ parm.traps = G_define_standard_option(G_OPT_R_INPUT);
+ parm.traps->key = "flow_control";
+ parm.traps->required = NO;
+ parm.traps->description =
+ _("Name of flow controls raster map (permeability ratio 0-1)");
+ parm.traps->guisection = _("Input");
+
+ parm.observation = G_define_standard_option(G_OPT_V_INPUT);
+ parm.observation->key = "observation";
+ parm.observation->required = NO;
+ parm.observation->label =
+ _("Name of sampling locations vector points map");
+ parm.observation->guisection = _("Input");
+
+ parm.depth = G_define_standard_option(G_OPT_R_OUTPUT);
+ parm.depth->key = "depth";
+ parm.depth->required = NO;
+ parm.depth->description = _("Name for output water depth raster map [m]");
+ parm.depth->guisection = _("Output");
+
+ parm.disch = G_define_standard_option(G_OPT_R_OUTPUT);
+ parm.disch->key = "discharge";
+ parm.disch->required = NO;
+ parm.disch->description = _("Name for output water discharge raster map [m3/s]");
+ parm.disch->guisection = _("Output");
+
+ parm.err = G_define_standard_option(G_OPT_R_OUTPUT);
+ parm.err->key = "error";
+ parm.err->required = NO;
+ parm.err->description = _("Name for output simulation error raster map [m]");
+ parm.err->guisection = _("Output");
+
+ parm.outwalk = G_define_standard_option(G_OPT_V_OUTPUT);
+ parm.outwalk->key = "walkers_output";
+ parm.outwalk->required = NO;
+ parm.outwalk->label =
+ _("Base name of the output walkers vector points map");
+ parm.outwalk->guisection = _("Output");
+
+ parm.logfile = G_define_standard_option(G_OPT_F_OUTPUT);
+ parm.logfile->key = "logfile";
+ parm.logfile->required = NO;
+ parm.logfile->description =
+ _("Name for sampling points output text file. For each observation vector point the time series of sediment transport is stored.");
+ parm.logfile->guisection = _("Output");
+
+ parm.nwalk = G_define_option();
+ parm.nwalk->key = "nwalkers";
+ parm.nwalk->type = TYPE_INTEGER;
+ parm.nwalk->required = NO;
+ parm.nwalk->description =
+ _("Number of walkers, default is twice the number of cells");
+ parm.nwalk->guisection = _("Parameters");
+
+ parm.niter = G_define_option();
+ parm.niter->key = "niterations";
+ parm.niter->type = TYPE_INTEGER;
+ parm.niter->answer = NITER;
+ parm.niter->required = NO;
+ parm.niter->description = _("Time used for iterations [minutes]");
+ parm.niter->guisection = _("Parameters");
+
+ parm.outiter = G_define_option();
+ parm.outiter->key = "output_step";
+ parm.outiter->type = TYPE_INTEGER;
+ parm.outiter->answer = ITEROUT;
+ parm.outiter->required = NO;
+ parm.outiter->description =
+ _("Time interval for creating output maps [minutes]");
+ parm.outiter->guisection = _("Parameters");
+
+/*
+ parm.density = G_define_option();
+ parm.density->key = "density";
+ parm.density->type = TYPE_INTEGER;
+ parm.density->answer = DENSITY;
+ parm.density->required = NO;
+ parm.density->description = _("Density of output walkers");
+ parm.density->guisection = _("Parameters");
+*/
+
+ parm.diffc = G_define_option();
+ parm.diffc->key = "diffusion_coeff";
+ parm.diffc->type = TYPE_DOUBLE;
+ parm.diffc->answer = DIFFC;
+ parm.diffc->required = NO;
+ parm.diffc->description = _("Water diffusion constant");
+ parm.diffc->guisection = _("Parameters");
+
+ parm.hmax = G_define_option();
+ parm.hmax->key = "hmax";
+ parm.hmax->type = TYPE_DOUBLE;
+ parm.hmax->answer = HMAX;
+ parm.hmax->required = NO;
+ parm.hmax->label =
+ _("Threshold water depth [m]");
+ parm.hmax->description = _("Diffusion increases after this water depth is reached");
+ parm.hmax->guisection = _("Parameters");
+
+ parm.halpha = G_define_option();
+ parm.halpha->key = "halpha";
+ parm.halpha->type = TYPE_DOUBLE;
+ parm.halpha->answer = HALPHA;
+ parm.halpha->required = NO;
+ parm.halpha->description = _("Diffusion increase constant");
+ parm.halpha->guisection = _("Parameters");
+
+ parm.hbeta = G_define_option();
+ parm.hbeta->key = "hbeta";
+ parm.hbeta->type = TYPE_DOUBLE;
+ parm.hbeta->answer = HBETA;
+ parm.hbeta->required = NO;
+ parm.hbeta->description =
+ _("Weighting factor for water flow velocity vector");
+ parm.hbeta->guisection = _("Parameters");
+
+#ifdef PARALLEL
+ parm.threads = G_define_option();
+ parm.threads->key = "threads";
+ parm.threads->type = TYPE_INTEGER;
+ parm.threads->answer = NUM_THREADS;
+ parm.threads->required = NO;
+ parm.threads->description =
+ _("Number of threads which will be used for parallel compute");
+ parm.threads->guisection = _("Parameters");
+#endif
+
+ flag.tserie = G_define_flag();
+ flag.tserie->key = 't';
+ flag.tserie->description = _("Time-series output");
+ flag.tserie->guisection = _("Output");
+
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+
+#ifdef PARALLEL
+ sscanf(parm.threads->answer, "%d", &threads);
+ if (threads < 1)
+ {
+ G_warning(_("<%d> is not valid number of threads. Number of threads will be set on <%d>"),
+ threads, abs(threads));
+ threads = abs(threads);
+ }
+ omp_set_num_threads(threads);
+ G_message(_("Number of treads <%d>"),threads);
+#endif
+
+ G_get_set_window(&cellhd);
+
+ WaterParams_init(&wp);
+
+ wp.conv = G_database_units_to_meters_factor();
+
+ G_debug(3, "Conversion factor is set to: %f", wp.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); */
+ 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", wp.xmax, wp.ymax);
+
+ wp.ts = flag.tserie->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", &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");
+
+ /* if no rain map input, then: */
+ if (parm.rain->answer == NULL) {
+ /*Check for Rain Unique Value Input */
+ /* if no rain unique value input */
+ if (parm.rainval->answer == NULL) {
+ /*No rain input so use default */
+ sscanf(RAINVAL, "%lf", &wp.rain_val);
+ /* if rain unique input exist, load it */
+ }
+ else {
+ /*Unique value input only */
+ 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) {
+ wp.rain_val = -999.99;
+ }
+ else {
+ /*both map and unique value exist */
+ /*Choose the map, discard the unique value */
+ wp.rain_val = -999.99;
+ }
+ }
+ /* Report the final value of 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) {
+ /*Check for Manin Unique Value Input */
+ /* if no Mannings unique value input */
+ if (parm.maninval->answer == NULL) {
+ /*No Mannings input so use default */
+ sscanf(MANINVAL, "%lf", &wp.manin_val);
+ /* if Mannings unique input value exists, load it */
+ }
+ else {
+ /*Unique value input only */
+ 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) {
+ wp.manin_val = -999.99;
+ }
+ else {
+ /*both map and unique value exist */
+ /*Choose map, discard the unique value */
+ wp.manin_val = -999.99;
+ }
+ }
+ /* Report the final value of 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) {
+ /*Check for Infil Unique Value Input */
+ /*if no infiltration unique value input */
+ if (parm.infilval->answer == NULL) {
+ /*No infiltration unique value so use default */
+ sscanf(INFILVAL, "%lf", &wp.infil_val);
+ /* if infiltration unique value exists, load it */
+ }
+ else {
+ /*unique value input only */
+ 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) {
+ wp.infil_val = -999.99;
+ }
+ else {
+ /*both map and unique value exist */
+ /*Choose map, discard the unique value */
+ wp.infil_val = -999.99;
+ }
+ }
+ /* Report the final value of 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 */
+ 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) {
+ 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", &wp.maxwa);
+ wp.rwalk = (double)wp.maxwa;
+ }
+
+ /* rwalk = (double) maxwa; */
+
+ if (wp.conv != 1.0)
+ G_message(_("Using metric conversion factor %f, step=%f"), wp.conv,
+ wp.step);
+
+ 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();
+
+ G_debug(1, "seeding randoms");
+ G_srand48(rand1);
+ grad_check();
+#ifdef PARALLEL
+ printTimeDiff("M2");
+#endif
+ main_loop();
+#ifdef PARALLEL
+ printTimeDiff("M3");
+#endif
+
+ /* Exit with Success */
+ exit(EXIT_SUCCESS);
+}
Added: grass-addons/grass7/raster/r.sim.water.mp/observation_points.c
===================================================================
--- grass-addons/grass7/raster/r.sim.water.mp/observation_points.c (rev 0)
+++ grass-addons/grass7/raster/r.sim.water.mp/observation_points.c 2017-04-28 04:43:41 UTC (rev 70973)
@@ -0,0 +1,138 @@
+/* obervation_points.c (simlib), 14.mar.2011, SG */
+
+#include <grass/glocale.h>
+#include <grass/vector.h>
+#include "waterglobs.h"
+
+
+static void init_points(struct _points *, int);
+static void realloc_points(struct _points *, int);
+static void insert_next_point(struct _points *p, double x, double y, int cat);
+
+/* ************************************************************************* */
+
+void create_observation_points()
+{
+ int if_log = 0;
+ struct Map_info Map;
+ struct line_pnts *pts;
+ struct line_cats *cts;
+ double x, y;
+ int type, cat, i;
+ struct Cell_head cellhd;
+
+ if(observation != NULL)
+ if_log += 1;
+
+ 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, observation, "") < 0)
+ G_fatal_error(_("Unable to open vector map <%s>"), observation);
+
+ Vect_rewind(&Map);
+
+ pts = Vect_new_line_struct();
+ cts = Vect_new_cats_struct();
+
+ /* 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);
+
+ if(type == -2) {
+ break;
+ }
+
+ if(type == -1) {
+ Vect_close(&Map);
+ G_fatal_error(_("Unable to read points from map %s"), observation);
+ }
+
+ if(type == GV_POINT) {
+ x = pts->x[0];
+ y = pts->y[0];
+ cat = cts->cat[0];
+
+ /* Check region bounds befor inserting point */
+ if(x <= cellhd.east && x >= cellhd.west && y <= cellhd.north && y >= cellhd.south) {
+ insert_next_point(&points, x, y, cat);
+ }
+ }
+ }
+
+ Vect_close(&Map);
+
+ /* Open the logfile */
+ points.output = fopen(logfile, "w");
+
+ if(points.output == NULL)
+ G_fatal_error(_("Unable to open observation logfile %s for writing"), logfile);
+
+ points.is_open = 1;
+
+ fprintf(points.output, "STEP ");
+
+ /* Write the vector cats as header information in the logfile */
+ for(i = 0; i < points.npoints; i++)
+ {
+ fprintf(points.output, "CAT%.4d ", points.cats[i]);
+ }
+ fprintf(points.output, "\n");
+
+
+ return;
+}
+
+/* ************************************************************************* */
+
+void init_points(struct _points *p, int size)
+{
+ p->x = (double*)G_calloc(size, sizeof(double));
+ p->y = (double*)G_calloc(size, sizeof(double));
+ p->cats = (int*)G_calloc(size, sizeof(int));
+ p->npoints = 0;
+ p->npoints_alloc = size;
+ p->output = NULL;
+ p->is_open = 0;
+}
+
+/* ************************************************************************* */
+
+void realloc_points(struct _points *p, int add_size)
+{
+ p->x = (double*)G_realloc(p->x, (p->npoints_alloc + add_size) * sizeof(double));
+ p->y = (double*)G_realloc(p->y, (p->npoints_alloc + add_size) * sizeof(double));
+ p->cats = (int*)G_realloc(p->cats, (p->npoints_alloc + add_size) * sizeof(int));
+ p->npoints_alloc += add_size;
+}
+
+/* ************************************************************************* */
+
+void insert_next_point(struct _points *p, double x, double y, int cat)
+{
+ if(p->npoints == p->npoints_alloc)
+ realloc_points(p, 128);
+
+ G_debug(3, "Insert point %g %g %i id %i\n", x, y, cat, p->npoints);
+
+ p->x[p->npoints] = x;
+ p->y[p->npoints] = y;
+ p->cats[p->npoints] = cat;
+ p->npoints++;
+}
Added: grass-addons/grass7/raster/r.sim.water.mp/output.c
===================================================================
--- grass-addons/grass7/raster/r.sim.water.mp/output.c (rev 0)
+++ grass-addons/grass7/raster/r.sim.water.mp/output.c 2017-04-28 04:43:41 UTC (rev 70973)
@@ -0,0 +1,764 @@
+/* output.c (simlib), 20.nov.2002, JH */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/bitmap.h>
+#include <grass/linkm.h>
+#include <grass/vector.h>
+#include <grass/glocale.h>
+
+#include "waterglobs.h"
+
+
+static void output_walker_as_vector(int tt_minutes, int ndigit, struct TimeStamp *timestamp);
+
+/* This function was added by Soeren 8. Mar 2011 */
+/* It replaces the site walker output implementation */
+/* Only the 3d coordinates of the walker are stored. */
+void output_walker_as_vector(int tt_minutes, int ndigit, struct TimeStamp *timestamp)
+{
+ char buf[GNAME_MAX + 10];
+ char *outwalk_time = NULL;
+ double x, y, z;
+ struct Map_info Out;
+ struct line_pnts *Points;
+ struct line_cats *Cats;
+ int i;
+
+ if (outwalk != NULL) {
+
+ /* In case of time series we extent the output name with the time value */
+ if (ts == 1) {
+ G_snprintf(buf, sizeof(buf), "%s_%.*d", outwalk, ndigit, tt_minutes);
+ outwalk_time = G_store(buf);
+ if (Vect_open_new(&Out, outwalk_time, WITH_Z) < 0)
+ G_fatal_error(_("Unable to create vector map <%s>"),
+ outwalk_time);
+ G_message("Writing %i walker into vector file %s", nstack, outwalk_time);
+ }
+ else {
+ if (Vect_open_new(&Out, outwalk, WITH_Z) < 0)
+ G_fatal_error(_("Unable to create vector map <%s>"), outwalk);
+ G_message("Writing %i walker into vector file %s", nstack, outwalk);
+ }
+
+ Points = Vect_new_line_struct();
+ Cats = Vect_new_cats_struct();
+
+ for (i = 0; i < nstack; i++) {
+ x = stack[i][0];
+ y = stack[i][1];
+ z = stack[i][2];
+
+ Vect_reset_line(Points);
+ Vect_reset_cats(Cats);
+
+ Vect_cat_set(Cats, 1, i + 1);
+ Vect_append_point(Points, x, y, z);
+ Vect_write_line(&Out, GV_POINT, Points, Cats);
+ }
+ /* Close vector file */
+ Vect_close(&Out);
+
+ Vect_destroy_line_struct(Points);
+ Vect_destroy_cats_struct(Cats);
+ if (ts == 1)
+ G_write_vector_timestamp(outwalk_time, "1", timestamp);
+ else
+ G_write_vector_timestamp(outwalk, "1", timestamp);
+ }
+
+ return;
+}
+
+/* Soeren 8. Mar 2011 TODO:
+ * This function needs to be refractured and splittet into smaller parts */
+int output_data(int tt, double ft)
+{
+
+ FCELL *depth_cell, *disch_cell, *err_cell;
+ FCELL *conc_cell, *flux_cell, *erdep_cell;
+ int depth_fd, disch_fd, err_fd;
+ int conc_fd, flux_fd, erdep_fd;
+ int i, iarc, j;
+ float gsmax = 0, dismax = 0., gmax = 0., ermax = -1.e+12, ermin = 1.e+12;
+ struct Colors colors;
+ struct History hist, hist1; /* hist2, hist3, hist4, hist5 */
+ struct TimeStamp timestamp;
+ char *depth0 = NULL, *disch0 = NULL, *err0 = NULL;
+ char *conc0 = NULL, *flux0 = NULL;
+ char *erdep0 = NULL;
+ const char *mapst = NULL;
+ char *type;
+ char buf[GNAME_MAX + 10];
+ char timestamp_buf[15];
+ int ndigit;
+ int timemin;
+ int tt_minutes;
+ FCELL dat1, dat2;
+ float a1, a2;
+
+ timemin = (int)(timesec / 60. + 0.5);
+ ndigit = 2;
+ /* more compact but harder to read:
+ ndigit = (int)floor(log10(timesec)) + 2 */
+ if (timemin >= 100)
+ ndigit = 3;
+ if (timemin >= 1000)
+ ndigit = 4;
+ if (timemin >= 10000)
+ ndigit = 5;
+
+ /* Convert to minutes */
+ tt_minutes = (int)(tt / 60. + 0.5);
+
+ /* Create timestamp */
+ sprintf(timestamp_buf, "%d minutes", tt_minutes);
+ G_scan_timestamp(×tamp, timestamp_buf);
+
+ /* Write the output walkers */
+ output_walker_as_vector(tt_minutes, ndigit, ×tamp);
+
+ /* 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,
+ Rast_window_rows());
+ if (mx != Rast_window_cols())
+ G_fatal_error("OOPS: cols changed from %d to %d", my,
+ Rast_window_cols());
+
+ if (depth) {
+ depth_cell = Rast_allocate_f_buf();
+ if (ts == 1) {
+ G_snprintf(buf, sizeof(buf), "%s.%.*d", depth, ndigit, tt_minutes);
+ depth0 = G_store(buf);
+ depth_fd = Rast_open_fp_new(depth0);
+ }
+ else
+ depth_fd = Rast_open_fp_new(depth);
+ }
+
+ if (disch) {
+ disch_cell = Rast_allocate_f_buf();
+ if (ts == 1) {
+ G_snprintf(buf, sizeof(buf),"%s.%.*d", disch, ndigit, tt_minutes);
+ disch0 = G_store(buf);
+ disch_fd = Rast_open_fp_new(disch0);
+ }
+ else
+ disch_fd = Rast_open_fp_new(disch);
+ }
+
+ if (err) {
+ err_cell = Rast_allocate_f_buf();
+ if (ts == 1) {
+ G_snprintf(buf, sizeof(buf), "%s.%.*d", err, ndigit, tt_minutes);
+ err0 = G_store(buf);
+ err_fd = Rast_open_fp_new(err0);
+ }
+ else
+ err_fd = Rast_open_fp_new(err);
+ }
+
+ if (conc) {
+ conc_cell = Rast_allocate_f_buf();
+ if (ts == 1) {
+ G_snprintf(buf, sizeof(buf), "%s.%.*d", conc, ndigit, tt_minutes);
+ conc0 = G_store(buf);
+ conc_fd = Rast_open_fp_new(conc0);
+ }
+ else
+ conc_fd = Rast_open_fp_new(conc);
+ }
+
+ if (flux) {
+ flux_cell = Rast_allocate_f_buf();
+ if (ts == 1) {
+ G_snprintf(buf, sizeof(buf), "%s.%.*d", flux, ndigit, tt_minutes);
+ flux0 = G_store(buf);
+ flux_fd = Rast_open_fp_new(flux0);
+ }
+ else
+ flux_fd = Rast_open_fp_new(flux);
+ }
+
+ if (erdep) {
+ erdep_cell = Rast_allocate_f_buf();
+ if (ts == 1) {
+ G_snprintf(buf, sizeof(buf), "%s.%.*d", erdep, ndigit, tt_minutes);
+ erdep0 = G_store(buf);
+ erdep_fd = Rast_open_fp_new(erdep0);
+ }
+ else
+ erdep_fd = Rast_open_fp_new(erdep);
+ }
+
+ for (iarc = 0; iarc < my; iarc++) {
+ i = my - iarc - 1;
+ if (depth) {
+ for (j = 0; j < mx; j++) {
+ if (zz[i][j] == UNDEF || gama[i][j] == UNDEF)
+ Rast_set_f_null_value(depth_cell + j, 1);
+ else {
+ a1 = pow(gama[i][j], 3. / 5.);
+ depth_cell[j] = (FCELL) a1; /* add conv? */
+ gmax = amax1(gmax, a1);
+ }
+ }
+ Rast_put_f_row(depth_fd, depth_cell);
+ }
+
+ if (disch) {
+ for (j = 0; j < mx; j++) {
+ if (zz[i][j] == UNDEF || gama[i][j] == UNDEF ||
+ cchez[i][j] == UNDEF)
+ Rast_set_f_null_value(disch_cell + j, 1);
+ else {
+ a2 = step * gama[i][j] * cchez[i][j]; /* cchez incl. sqrt(sinsl) */
+ disch_cell[j] = (FCELL) a2; /* add conv? */
+ dismax = amax1(dismax, a2);
+ }
+ }
+ Rast_put_f_row(disch_fd, disch_cell);
+ }
+
+ if (err) {
+ for (j = 0; j < mx; j++) {
+ if (zz[i][j] == UNDEF || gammas[i][j] == UNDEF)
+ Rast_set_f_null_value(err_cell + j, 1);
+ else {
+ err_cell[j] = (FCELL) gammas[i][j];
+ gsmax = amax1(gsmax, gammas[i][j]); /* add conv? */
+ }
+ }
+ Rast_put_f_row(err_fd, err_cell);
+ }
+
+ if (conc) {
+ for (j = 0; j < mx; j++) {
+ if (zz[i][j] == UNDEF || gama[i][j] == UNDEF)
+ Rast_set_f_null_value(conc_cell + j, 1);
+ else {
+ conc_cell[j] = (FCELL) gama[i][j];
+ /* gsmax = amax1(gsmax, gama[i][j]); */
+ }
+ }
+ Rast_put_f_row(conc_fd, conc_cell);
+ }
+
+ if (flux) {
+ for (j = 0; j < mx; j++) {
+ if (zz[i][j] == UNDEF || gama[i][j] == UNDEF ||
+ slope[i][j] == UNDEF)
+ Rast_set_f_null_value(flux_cell + j, 1);
+ else {
+ a2 = gama[i][j] * slope[i][j];
+ flux_cell[j] = (FCELL) a2;
+ dismax = amax1(dismax, a2);
+ }
+ }
+ Rast_put_f_row(flux_fd, flux_cell);
+ }
+
+ if (erdep) {
+ for (j = 0; j < mx; j++) {
+ if (zz[i][j] == UNDEF || er[i][j] == UNDEF)
+ Rast_set_f_null_value(erdep_cell + j, 1);
+ else {
+ erdep_cell[j] = (FCELL) er[i][j];
+ ermax = amax1(ermax, er[i][j]);
+ ermin = amin1(ermin, er[i][j]);
+ }
+ }
+ Rast_put_f_row(erdep_fd, erdep_cell);
+ }
+ }
+
+ if (depth)
+ Rast_close(depth_fd);
+ if (disch)
+ Rast_close(disch_fd);
+ if (err)
+ Rast_close(err_fd);
+ if (conc)
+ Rast_close(conc_fd);
+ if (flux)
+ Rast_close(flux_fd);
+ if (erdep)
+ Rast_close(erdep_fd);
+
+ if (depth) {
+
+ Rast_init_colors(&colors);
+
+ dat1 = (FCELL) 0.;
+ dat2 = (FCELL) 0.001;
+ Rast_add_f_color_rule(&dat1, 255, 255, 255, &dat2, 255, 255, 0,
+ &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) 0.05;
+ Rast_add_f_color_rule(&dat1, 255, 255, 0, &dat2, 0, 255, 255,
+ &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) 0.1;
+ Rast_add_f_color_rule(&dat1, 0, 255, 255, &dat2, 0, 127, 255,
+ &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) 0.5;
+ Rast_add_f_color_rule(&dat1, 0, 127, 255, &dat2, 0, 0, 255,
+ &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) gmax;
+ Rast_add_f_color_rule(&dat1, 0, 0, 255, &dat2, 0, 0, 0, &colors);
+
+
+ if (ts == 1) {
+ if ((mapst = G_find_file("fcell", depth0, "")) == NULL)
+ G_fatal_error(_("FP raster map <%s> not found"), depth0);
+ Rast_write_colors(depth0, mapst, &colors);
+ Rast_quantize_fp_map_range(depth0, mapst, 0., (FCELL) gmax, 0,
+ (CELL) gmax);
+ Rast_free_colors(&colors);
+ }
+ else {
+ if ((mapst = G_find_file("fcell", depth, "")) == NULL)
+ G_fatal_error(_("FP raster map <%s> not found"), depth);
+ Rast_write_colors(depth, mapst, &colors);
+ Rast_quantize_fp_map_range(depth, mapst, 0., (FCELL) gmax, 0,
+ (CELL) gmax);
+ Rast_free_colors(&colors);
+ }
+ }
+
+ if (disch) {
+
+ Rast_init_colors(&colors);
+
+ dat1 = (FCELL) 0.;
+ dat2 = (FCELL) 0.0005;
+ Rast_add_f_color_rule(&dat1, 255, 255, 255, &dat2, 255, 255, 0,
+ &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) 0.005;
+ Rast_add_f_color_rule(&dat1, 255, 255, 0, &dat2, 0, 255, 255,
+ &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) 0.05;
+ Rast_add_f_color_rule(&dat1, 0, 255, 255, &dat2, 0, 127, 255,
+ &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) 0.1;
+ Rast_add_f_color_rule(&dat1, 0, 127, 255, &dat2, 0, 0, 255,
+ &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) dismax;
+ Rast_add_f_color_rule(&dat1, 0, 0, 255, &dat2, 0, 0, 0, &colors);
+
+ if (ts == 1) {
+ if ((mapst = G_find_file("cell", disch0, "")) == NULL)
+ G_fatal_error(_("Raster map <%s> not found"), disch0);
+ Rast_write_colors(disch0, mapst, &colors);
+ Rast_quantize_fp_map_range(disch0, mapst, 0., (FCELL) dismax, 0,
+ (CELL) dismax);
+ Rast_free_colors(&colors);
+ }
+ else {
+
+ if ((mapst = G_find_file("cell", disch, "")) == NULL)
+ G_fatal_error(_("Raster map <%s> not found"), disch);
+ Rast_write_colors(disch, mapst, &colors);
+ Rast_quantize_fp_map_range(disch, mapst, 0., (FCELL) dismax, 0,
+ (CELL) dismax);
+ Rast_free_colors(&colors);
+ }
+ }
+
+ if (flux) {
+
+ Rast_init_colors(&colors);
+
+ dat1 = (FCELL) 0.;
+ dat2 = (FCELL) 0.001;
+ Rast_add_f_color_rule(&dat1, 255, 255, 255, &dat2, 255, 255, 0,
+ &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) 0.1;
+ Rast_add_f_color_rule(&dat1, 255, 255, 0, &dat2, 255, 127, 0,
+ &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) 1.;
+ Rast_add_f_color_rule(&dat1, 255, 127, 0, &dat2, 191, 127, 63,
+ &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) dismax;
+ Rast_add_f_color_rule(&dat1, 191, 127, 63, &dat2, 0, 0, 0,
+ &colors);
+
+ if (ts == 1) {
+ if ((mapst = G_find_file("cell", flux0, "")) == NULL)
+ G_fatal_error(_("Raster map <%s> not found"), flux0);
+ Rast_write_colors(flux0, mapst, &colors);
+ Rast_quantize_fp_map_range(flux0, mapst, 0., (FCELL) dismax, 0,
+ (CELL) dismax);
+ Rast_free_colors(&colors);
+ }
+ else {
+
+ if ((mapst = G_find_file("cell", flux, "")) == NULL)
+ G_fatal_error(_("Raster map <%s> not found"), flux);
+ Rast_write_colors(flux, mapst, &colors);
+ Rast_quantize_fp_map_range(flux, mapst, 0., (FCELL) dismax, 0,
+ (CELL) dismax);
+ Rast_free_colors(&colors);
+ }
+ }
+
+ if (erdep) {
+
+ Rast_init_colors(&colors);
+
+ dat1 = (FCELL) ermax;
+ dat2 = (FCELL) 0.1;
+ Rast_add_f_color_rule(&dat1, 0, 0, 0, &dat2, 0, 0, 255, &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) 0.01;
+ Rast_add_f_color_rule(&dat1, 0, 0, 255, &dat2, 0, 191, 191,
+ &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) 0.0001;
+ Rast_add_f_color_rule(&dat1, 0, 191, 191, &dat2, 170, 255, 255,
+ &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) 0.;
+ Rast_add_f_color_rule(&dat1, 170, 255, 255, &dat2, 255, 255, 255,
+ &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) - 0.0001;
+ Rast_add_f_color_rule(&dat1, 255, 255, 255, &dat2, 255, 255, 0,
+ &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) - 0.01;
+ Rast_add_f_color_rule(&dat1, 255, 255, 0, &dat2, 255, 127, 0,
+ &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) - 0.1;
+ Rast_add_f_color_rule(&dat1, 255, 127, 0, &dat2, 255, 0, 0,
+ &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) ermin;
+ Rast_add_f_color_rule(&dat1, 255, 0, 0, &dat2, 255, 0, 255,
+ &colors);
+
+ if (ts == 1) {
+ if ((mapst = G_find_file("cell", erdep0, "")) == NULL)
+ G_fatal_error(_("Raster map <%s> not found"), erdep0);
+ Rast_write_colors(erdep0, mapst, &colors);
+ Rast_quantize_fp_map_range(erdep0, mapst, (FCELL) ermin,
+ (FCELL) ermax, (CELL) ermin,
+ (CELL) ermax);
+ Rast_free_colors(&colors);
+
+ type = "raster";
+ Rast_short_history(erdep0, type, &hist1);
+ Rast_append_format_history(
+ &hist1, "The sediment flux file is %s", flux0);
+ Rast_command_history(&hist1);
+ Rast_write_history(erdep0, &hist1);
+ }
+ else {
+
+ if ((mapst = G_find_file("cell", erdep, "")) == NULL)
+ G_fatal_error(_("Raster map <%s> not found"), erdep);
+ Rast_write_colors(erdep, mapst, &colors);
+ Rast_quantize_fp_map_range(erdep, mapst, (FCELL) ermin,
+ (FCELL) ermax, (CELL) ermin,
+ (CELL) ermax);
+ Rast_free_colors(&colors);
+
+ type = "raster";
+ Rast_short_history(erdep, type, &hist1);
+ Rast_append_format_history(
+ &hist1, "The sediment flux file is %s", flux);
+ Rast_command_history(&hist1);
+ Rast_write_history(erdep, &hist1);
+ }
+ }
+
+ /* history section */
+ if (depth) {
+ type = "raster";
+ if (ts == 0) {
+ mapst = G_find_file("cell", depth, "");
+ if (mapst == NULL) {
+ G_warning(_("Raster map <%s> not found"), depth);
+ return -1;
+ }
+ Rast_short_history(depth, type, &hist);
+ }
+ else
+ Rast_short_history(depth0, type, &hist);
+
+ Rast_append_format_history(
+ &hist, "init.walk=%d, maxwalk=%d, remaining walkers=%d",
+ nwalk, maxwa, nwalka);
+ Rast_append_format_history(
+ &hist, "duration (sec.)=%d, time-serie iteration=%d",
+ timesec, tt);
+ Rast_append_format_history(
+ &hist, "written deltap=%f, mean vel.=%f",
+ deltap, vmean);
+ Rast_append_format_history(
+ &hist, "mean source (si)=%e, mean infil=%e",
+ si0, infmean);
+
+ Rast_format_history(&hist, HIST_DATSRC_1, "input files: %s %s %s",
+ elevin, dxin, dyin);
+ Rast_format_history(&hist, HIST_DATSRC_2, "input files: %s %s %s",
+ rain, infil, manin);
+
+ Rast_command_history(&hist);
+
+ if (ts == 1)
+ Rast_write_history(depth0, &hist);
+ else
+ Rast_write_history(depth, &hist);
+
+ if (ts == 1)
+ G_write_raster_timestamp(depth0, ×tamp);
+ else
+ G_write_raster_timestamp(depth, ×tamp);
+ }
+
+ if (disch) {
+ type = "raster";
+ if (ts == 0) {
+ mapst = G_find_file("cell", disch, "");
+ if (mapst == NULL)
+ G_fatal_error(_("Raster map <%s> not found"), disch);
+ Rast_short_history(disch, type, &hist);
+ }
+ else
+ Rast_short_history(disch0, type, &hist);
+
+ Rast_append_format_history(
+ &hist,"init.walkers=%d, maxwalk=%d, rem. walkers=%d",
+ nwalk, maxwa, nwalka);
+ Rast_append_format_history(
+ &hist, "duration (sec.)=%d, time-serie iteration=%d",
+ timesec, tt);
+ Rast_append_format_history(
+ &hist, "written deltap=%f, mean vel.=%f",
+ deltap, vmean);
+ Rast_append_format_history(
+ &hist, "mean source (si)=%e, mean infil=%e",
+ si0, infmean);
+
+ Rast_format_history(&hist, HIST_DATSRC_1, "input files: %s %s %s",
+ elevin, dxin, dyin);
+ Rast_format_history(&hist, HIST_DATSRC_2, "input files: %s %s %s",
+ rain, infil, manin);
+
+ Rast_command_history(&hist);
+
+ if (ts == 1)
+ Rast_write_history(disch0, &hist);
+ else
+ Rast_write_history(disch, &hist);
+
+ if (ts == 1)
+ G_write_raster_timestamp(disch0, ×tamp);
+ else
+ G_write_raster_timestamp(disch, ×tamp);
+ }
+
+ if (flux) {
+ type = "raster";
+ if (ts == 0) {
+ mapst = G_find_file("cell", flux, "");
+ if (mapst == NULL)
+ G_fatal_error(_("Raster map <%s> not found"), flux);
+ Rast_short_history(flux, type, &hist);
+ }
+ else
+ Rast_short_history(flux0, type, &hist);
+
+ Rast_append_format_history(
+ &hist, "init.walk=%d, maxwalk=%d, remaining walkers=%d",
+ nwalk, maxwa, nwalka);
+ Rast_append_format_history(
+ &hist, "duration (sec.)=%d, time-serie iteration=%d",
+ timesec, tt);
+ Rast_append_format_history(
+ &hist, "written deltap=%f, mean vel.=%f",
+ deltap, vmean);
+ Rast_append_format_history(
+ &hist, "mean source (si)=%f", si0);
+
+ Rast_format_history(&hist, HIST_DATSRC_1, "input files: %s %s %s",
+ wdepth, dxin, dyin);
+ Rast_format_history(&hist, HIST_DATSRC_2, "input files: %s %s %s %s",
+ manin, detin, tranin, tauin);
+
+ Rast_command_history(&hist);
+
+ if (ts == 1)
+ Rast_write_history(flux0, &hist);
+ else
+ Rast_write_history(flux, &hist);
+
+ if (ts == 1)
+ G_write_raster_timestamp(flux0, ×tamp);
+ else
+ G_write_raster_timestamp(flux, ×tamp);
+ }
+
+ return 1;
+}
+
+
+int output_et()
+{
+
+ FCELL *tc_cell, *et_cell;
+ int tc_fd, et_fd;
+ int i, iarc, j;
+ float etmax = -1.e+12, etmin = 1.e+12;
+ float trc;
+ struct Colors colors;
+ const char *mapst = NULL;
+
+ /* char buf[GNAME_MAX + 10]; */
+ FCELL dat1, dat2;
+
+ /* float a1,a2; */
+
+ /* we write in the same region as we used for reading */
+
+ if (et) {
+ et_cell = Rast_allocate_f_buf();
+ /* if (ts == 1) {
+ sprintf(buf, "%s.%.*d", et, ndigit, tt);
+ et0 = G_store(buf);
+ et_fd = Rast_open_fp_new(et0);
+ }
+ else */
+ et_fd = Rast_open_fp_new(et);
+ }
+
+ if (tc) {
+ tc_cell = Rast_allocate_f_buf();
+ /* if (ts == 1) {
+ sprintf(buf, "%s.%.*d", tc, ndigit, tt);
+ tc0 = G_store(buf);
+ tc_fd = Rast_open_fp_new(tc0);
+ }
+ else */
+ tc_fd = Rast_open_fp_new(tc);
+ }
+
+ if (my != Rast_window_rows())
+ G_fatal_error("OOPS: rows changed from %d to %d", mx,
+ Rast_window_rows());
+ if (mx != Rast_window_cols())
+ G_fatal_error("OOPS: cols changed from %d to %d", my,
+ Rast_window_cols());
+
+ for (iarc = 0; iarc < my; iarc++) {
+ i = my - iarc - 1;
+ if (et) {
+ for (j = 0; j < mx; j++) {
+ if (zz[i][j] == UNDEF || er[i][j] == UNDEF)
+ Rast_set_f_null_value(et_cell + j, 1);
+ else {
+ et_cell[j] = (FCELL) er[i][j]; /* add conv? */
+ etmax = amax1(etmax, er[i][j]);
+ etmin = amin1(etmin, er[i][j]);
+ }
+ }
+ Rast_put_f_row(et_fd, et_cell);
+ }
+
+ if (tc) {
+ for (j = 0; j < mx; j++) {
+ if (zz[i][j] == UNDEF || sigma[i][j] == UNDEF ||
+ si[i][j] == UNDEF)
+ Rast_set_f_null_value(tc_cell + j, 1);
+ else {
+ if (sigma[i][j] == 0.)
+ trc = 0.;
+ else
+ trc = si[i][j] / sigma[i][j];
+ tc_cell[j] = (FCELL) trc;
+ /* gsmax = amax1(gsmax, trc); */
+ }
+ }
+ Rast_put_f_row(tc_fd, tc_cell);
+ }
+ }
+
+ if (tc)
+ Rast_close(tc_fd);
+
+ if (et)
+ Rast_close(et_fd);
+
+ if (et) {
+
+ Rast_init_colors(&colors);
+
+ dat1 = (FCELL) etmax;
+ dat2 = (FCELL) 0.1;
+ Rast_add_f_color_rule(&dat1, 0, 0, 0, &dat2, 0, 0, 255, &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) 0.01;
+ Rast_add_f_color_rule(&dat1, 0, 0, 255, &dat2, 0, 191, 191,
+ &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) 0.0001;
+ Rast_add_f_color_rule(&dat1, 0, 191, 191, &dat2, 170, 255, 255,
+ &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) 0.;
+ Rast_add_f_color_rule(&dat1, 170, 255, 255, &dat2, 255, 255, 255,
+ &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) - 0.0001;
+ Rast_add_f_color_rule(&dat1, 255, 255, 255, &dat2, 255, 255, 0,
+ &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) - 0.01;
+ Rast_add_f_color_rule(&dat1, 255, 255, 0, &dat2, 255, 127, 0,
+ &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) - 0.1;
+ Rast_add_f_color_rule(&dat1, 255, 127, 0, &dat2, 255, 0, 0,
+ &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) etmin;
+ Rast_add_f_color_rule(&dat1, 255, 0, 0, &dat2, 255, 0, 255,
+ &colors);
+
+ /* if (ts == 1) {
+ if ((mapst = G_find_file("cell", et0, "")) == NULL)
+ G_fatal_error(_("Raster map <%s> not found"), et0);
+ Rast_write_colors(et0, mapst, &colors);
+ Rast_quantize_fp_map_range(et0, mapst, (FCELL)etmin, (FCELL)etmax, (CELL)etmin, (CELL)etmax);
+ Rast_free_colors(&colors);
+ }
+ else { */
+
+ if ((mapst = G_find_file("cell", et, "")) == NULL)
+ G_fatal_error(_("Raster map <%s> not found"), et);
+ Rast_write_colors(et, mapst, &colors);
+ Rast_quantize_fp_map_range(et, mapst, (FCELL) etmin, (FCELL) etmax,
+ (CELL) etmin, (CELL) etmax);
+ Rast_free_colors(&colors);
+ /* } */
+ }
+
+ return 1;
+}
Added: grass-addons/grass7/raster/r.sim.water.mp/r.sim.water.mp.html
===================================================================
--- grass-addons/grass7/raster/r.sim.water.mp/r.sim.water.mp.html (rev 0)
+++ grass-addons/grass7/raster/r.sim.water.mp/r.sim.water.mp.html 2017-04-28 04:43:41 UTC (rev 70973)
@@ -0,0 +1,240 @@
+<h2>DESCRIPTION</h2>
+
+<i>r.sim.water</i> is a landscape scale simulation model
+of overland flow designed for spatially variable terrain, soil, cover
+and rainfall excess conditions. A 2D shallow water flow is described by
+the bivariate form of Saint Venant equations. The numerical solution is based
+on the concept of duality between the field and particle representation of
+the modeled quantity. Green's function Monte Carlo method, used to solve the equation,
+provides robustness necessary for spatially variable conditions and high
+resolutions (Mitas and Mitasova 1998). The key inputs of the model include
+elevation (<i>elevation</i> raster map), flow gradient vector given by
+first-order partial derivatives of elevation field (<i>dx</i> and <i>dy</i>
+raster maps), rainfall excess rate (<i>rain</i> raster map or <i>rain_value</i> single
+value) and a surface roughness coefficient given by Manning's n
+(<i>man</i> raster map or <i>man_value</i> single value). Partial
+derivatives raster maps can be computed along with interpolation of a DEM using
+the -d option in <a href="v.surf.rst.html">v.surf.rst</a> module. If elevation raster
+map is already provided, partial derivatives can be computed using
+<a href="r.slope.aspect.html">r.slope.aspect</a> module. Partial derivatives are used
+to determine the direction and magnitude of water flow velocity. To include a
+predefined direction of flow, map algebra can be used to replace terrain-derived
+partial derivatives with pre-defined partial derivatives in selected grid cells such
+as man-made channels, ditches or culverts. Equations (2) and (3) from
+<a href="http://www4.ncsu.edu/~hmitaso/gmslab/reports/cerl99/rep99.html">this report</a>
+can be used to compute partial derivates of the predefined flow using its direction given
+by aspect and slope.
+
+<p>
+The module automatically converts horizontal distances from feet to metric system using
+database/projection information. Rainfall excess is defined as rainfall intensity
+- infiltration rate and should be provided in [mm/hr].
+<!-- and can be computed using several available infiltration
+ models (e.g. Green-Ampt, Holtan, etc.). (<font color="#ff0000"> find
+ infiltration module in GRASS - topmodel, casc2d</font> )
+-->
+Rainfall intensities are usually available from meteorological stations.
+Infiltration rate depends on soil properties and land cover. It varies in space and time.
+For saturated soil and steady-state water flow it can be estimated using
+saturated hydraulic conductivity rates based on field measurements or using
+reference values which can be found in literature.
+Optionally, user can provide an overland flow infiltration rate map
+<i>infil</i> or a single value <i>infil_value</i> in [mm/hr] that control the rate of
+infiltration for the already flowing water, effectively reducing the flow depth and
+discharge.
+Overland flow can be further controlled by permeable check dams or similar type of structures,
+the user can provide a map of these structures and their permeability ratio
+in the map <i>flow_control</i> that defines the probability of particles to pass
+through the structure (the values will be 0-1).
+
+<p>
+Output includes a water depth raster map <i>depth</i> in [m], and a water discharge
+raster map <i>discharge</i> in [m3/s]. Error of the numerical solution can be analyzed using
+the <i>error</i> raster map (the resulting water depth is an average, and err is its RMSE).
+The output vector points map <i>output_walkers</i> can be used to analyze and visualize
+spatial distribution of walkers at different simulation times (note that
+the resulting water depth is based on the density of these walkers).
+<!--Number of the output walkers is controlled by the <i>density</i> parameter, which controls
+how many walkers used in simulation should be written into the output. -->
+<!-- from
+http://www.ing.unitn.it/~grass/conferences/GRASS2002/proceedings/proceedings/pdfs/Mitasova_Helena_2.pdf
+-->
+The spatial distribution of numerical error associated with path sampling solution can be
+analysed using the output error raster file [m]. This error is a function of the number
+of particles used in the simulation and can be reduced by increasing the number of walkers
+given by parameter <i>nwalkers</i>.
+<!--(<font color="#ff0000"> toto treba upresnit/zmenit, lebo nwalk ide prec</font>). -->
+Duration of simulation is controlled by the <i>niterations</i> parameter. The default value
+is 10 minutes, reaching the steady-state may require much longer time,
+depending on the time step, complexity of terrain, land cover and size of the area.
+Output walker, water depth and discharge maps can be saved during simulation using
+the time series flag <i>-t</i> and <i>output_step</i> parameter
+defining the time step in minutes for writing output files.
+Files are saved with a suffix representing time since the start of simulation in minutes
+(e.g. wdepth.05, wdepth.10).
+<!-- TODO add example of registering outputs into temporal framework when #2146 is fixed -->
+Monitoring of water depth at specific points is supported. A vector map with observation points and
+a path to a logfile must be provided. For each point in the vector map which is located in
+the computational region the water depth is logged each time step in the logfile. The logfile is
+organized as a table. A single header identifies the category number of the logged vector points.
+In case of invalid water depth data the value -1 is used.
+
+<p>
+Overland flow is routed based on partial derivatives of elevation
+field or other landscape features influencing water flow. Simulation
+equations include a diffusion term (<i>diffusion_coeff</i> parameter) which enables
+water flow to overcome elevation depressions or obstacles when water depth exceeds
+a threshold water depth value (<i>hmax)</i>, given in [m]. When it is reached,
+diffusion term increases as given by <i>halpha</i> and advection term
+(direction of flow) is given as "prevailing" direction of flow computed
+as average of flow directions from the previous <i>hbeta</i> number of grid cells.
+
+<h2>NOTES</h2>
+
+A 2D shallow water flow is described by the bivariate form of Saint
+Venant equations (e.g., Julien et al., 1995). The continuity of water
+flow relation is coupled with the momentum conservation equation and
+for a shallow water overland flow, the hydraulic radius is approximated
+by the normal flow depth. The system of equations is closed using the
+Manning's relation. Model assumes that the flow is close to the kinematic
+wave approximation, but we include a diffusion-like term to incorporate the
+impact of diffusive wave effects. Such an incorporation of diffusion
+in the water flow simulation is not new and a similar term has been obtained
+in derivations of diffusion-advection equations for overland flow, e.g.,
+by Lettenmeier and Wood, (1992). In our reformulation, we simplify the
+diffusion coefficient to a constant and we use a modified diffusion term.
+The diffusion constant which we have used is rather small (approximately
+one order of magnitude smaller than the reciprocal Manning's coefficient)
+and therefore the resulting flow is close to the kinematic regime. However,
+the diffusion term improves the kinematic solution, by overcoming small
+shallow pits common in digital elevation models (DEM) and by smoothing out
+the flow over slope discontinuities or abrupt changes in Manning's coefficient
+(e.g., due to a road, or other anthropogenic changes in elevations or cover).
+
+<p>
+<b>Green's function stochastic method of solution.</b><br>
+The Saint Venant equations are solved by a stochastic method called Monte Carlo
+(very similar to Monte Carlo methods in computational fluid dynamics or to
+quantum Monte Carlo approaches for solving the Schrodinger equation (Schmidt
+and Ceperley, 1992, Hammond et al., 1994; Mitas, 1996)). It is assumed
+that these equations are a representation of stochastic processes with
+diffusion and drift components (Fokker-Planck equations).
+
+<p>
+The Monte Carlo technique has several unique advantages which are
+becoming even more important due to new developments in computer technology.
+Perhaps one of the most significant Monte Carlo properties is robustness
+which enables us to solve the equations for complex cases, such as discontinuities
+in the coefficients of differential operators (in our case, abrupt slope
+or cover changes, etc). Also, rough solutions can be estimated rather
+quickly, which allows us to carry out preliminary quantitative studies
+or to rapidly extract qualitative trends by parameter scans. In addition,
+the stochastic methods are tailored to the new generation of computers
+as they provide scalability from a single workstation to large parallel
+machines due to the independence of sampling points. Therefore, the methods
+are useful both for everyday exploratory work using a desktop computer and
+for large, cutting-edge applications using high performance computing.
+
+<h2>EXAMPLE</h2>
+
+Spearfish region:
+
+<div class="code"><pre>
+g.region raster=elevation.10m -p
+r.slope.aspect elevation=elevation.10m dx=elev_dx dy=elev_dy
+
+# synthetic maps
+r.mapcalc "rain = if(elevation.10m, 5.0, null())"
+r.mapcalc "manning = if(elevation.10m, 0.05, null())"
+r.mapcalc "infilt = if(elevation.10m, 0.0, null())"
+
+# simulate
+r.sim.water elevation=elevation.10m dx=elev_dx dy=elev_dy \
+ rain=rain man=manning infil=infilt \
+ nwalkers=5000000 depth=depth
+</pre></div>
+
+<p>
+<center>
+<img src="r_sim_water.png" alt="r.sim.water generated depth map"><br>
+<i>Water depth map in the Spearfish (SD) area</i>
+</center>
+
+
+<h2>ERROR MESSAGES</h2>
+
+If the module fails with
+
+<div class="code"><pre>
+ERROR: nwalk (7000001) > maxw (7000000)!
+</pre></div>
+
+then a lower <em>nwalkers</em> parameter value has to be selected.
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="v.surf.rst.html">v.surf.rst</a>,
+<a href="r.slope.aspect.html">r.slope.aspect</a>,
+<a href="r.sim.sediment.html">r.sim.sediment</a>
+</em>
+
+<h2>AUTHORS</h2>
+
+Helena Mitasova, Lubos Mitas<br>
+North Carolina State University<br>
+<i><a href="mailto:hmitaso at unity.ncsu.edu">hmitaso at unity.ncsu.edu</a></i>
+
+<p>
+Jaroslav Hofierka<br>
+GeoModel, s.r.o. Bratislava, Slovakia<br>
+<i><a href="mailto:hofi at geomodel.sk">hofierka at geomodel.sk</a></i>
+
+<p>
+Chris Thaxton<br>
+North Carolina State University<br>
+<i><a href="mailto:csthaxto at unity.ncsu.edu">csthaxto at unity.ncsu.edu</a></i>
+
+<h2>REFERENCES</h2>
+
+<ul>
+<li> Mitasova, H., Thaxton, C., Hofierka, J., McLaughlin, R., Moore, A., Mitas L., 2004,
+<a href="http://www4.ncsu.edu/~hmitaso/gmslab/papers/II.6.8_Mitasova_044.pdf">
+Path sampling method for modeling overland water flow, sediment transport
+and short term terrain evolution in Open Source GIS.</a>
+In: C.T. Miller, M.W. Farthing, V.G. Gray, G.F. Pinder eds.,
+Proceedings of the XVth International Conference on Computational Methods in Water
+Resources (CMWR XV), June 13-17 2004, Chapel Hill, NC, USA, Elsevier, pp. 1479-1490.
+
+<li> Mitasova H, Mitas, L., 2000,
+<a href="http://www4.ncsu.edu/~hmitaso/gmslab/gisc00/duality.html">Modeling spatial
+processes in multiscale framework: exploring duality between particles and fields,</a>
+plenary talk at GIScience2000 conference, Savannah, GA.
+
+<li> Mitas, L., and Mitasova, H., 1998, Distributed soil erosion simulation
+for effective erosion prevention. Water Resources Research, 34(3), 505-516.
+
+<li> Mitasova, H., Mitas, L., 2001,
+<a href="http://www4.ncsu.edu/~hmitaso/gmslab/papers/LLEmiterev1.pdf">
+Multiscale soil erosion simulations for land use management,</a>
+In: Landscape erosion and landscape evolution modeling, Harmon R. and Doe W. eds.,
+Kluwer Academic/Plenum Publishers, pp. 321-347.
+
+<li> Hofierka, J, Mitasova, H., Mitas, L., 2002. GRASS and modeling landscape processes
+using duality between particles and fields. Proceedings of the Open source GIS -
+GRASS users conference 2002 - Trento, Italy, 11-13 September 2002.
+<a href="http://www.ing.unitn.it/~grass/conferences/GRASS2002/proceedings/proceedings/pdfs/Mitasova_Helena_2.pdf">PDF</a>
+
+<li> Hofierka, J., Knutova, M., 2015,
+Simulating aspects of a flash flood using the Monte Carlo method and
+GRASS GIS: a case study of the Malá Svinka Basin (Slovakia),
+Open Geosciences. Volume 7, Issue 1, ISSN (Online) 2391-5447, DOI:
+<a href="http://dx.doi.org/10.1515/geo-2015-0013">10.1515/geo-2015-0013</a>,
+April 2015
+
+<li> Neteler, M. and Mitasova, H., 2008,
+<a href="http://www.grassbook.org">Open Source GIS: A GRASS GIS Approach. Third Edition.</a>
+The International Series in Engineering and Computer Science: Volume 773. Springer New York Inc, p. 406.
+</ul>
+
+<p><i>Last changed: $Date: 2016-03-08 09:06:33 +0100 (Út, 08 bře 2016) $</i>
Added: grass-addons/grass7/raster/r.sim.water.mp/r_sim_water.png
===================================================================
(Binary files differ)
Property changes on: grass-addons/grass7/raster/r.sim.water.mp/r_sim_water.png
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: grass-addons/grass7/raster/r.sim.water.mp/random.c
===================================================================
--- grass-addons/grass7/raster/r.sim.water.mp/random.c (rev 0)
+++ grass-addons/grass7/raster/r.sim.water.mp/random.c 2017-04-28 04:43:41 UTC (rev 70973)
@@ -0,0 +1,61 @@
+/* random.c (simlib), 20.nov.2002, JH */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <grass/gis.h>
+#include <grass/bitmap.h>
+#include <grass/linkm.h>
+
+#include "waterglobs.h"
+
+double simwe_rand(void)
+{
+ return G_drand48();
+} /* ulec */
+
+
+double gasdev(void)
+{
+ /* Initialized data */
+
+ static int iset = 0;
+ static double gset = .1;
+
+ /* System generated locals */
+ double ret_val;
+
+ /* Local variables */
+ double r = 0., vv1, vv2, fac;
+
+ if (iset == 0) {
+ while (r >= 1. || r == 0.) {
+ vv1 = simwe_rand() * 2. - 1.;
+ vv2 = simwe_rand() * 2. - 1.;
+ r = vv1 * vv1 + vv2 * vv2;
+ }
+ fac = sqrt(log(r) * -2. / r);
+ gset = vv1 * fac;
+ ret_val = vv2 * fac;
+ iset = 1;
+ }
+ else {
+ ret_val = gset;
+ iset = 0;
+ }
+ return ret_val;
+} /* gasdev */
+
+void gasdev_for_paralel(double *x, double *y )
+{
+ double r = 0., vv1, vv2, fac;
+
+ while (r >= 1. || r == 0.) {
+ vv1 = simwe_rand() * 2. - 1.;
+ vv2 = simwe_rand() * 2. - 1.;
+ r = vv1 * vv1 + vv2 * vv2;
+ }
+ fac = sqrt(log(r) * -2. / r);
+ (*y) = vv1 * fac;
+ (*x) = vv2 * fac;
+}
Added: grass-addons/grass7/raster/r.sim.water.mp/simlib.h
===================================================================
--- grass-addons/grass7/raster/r.sim.water.mp/simlib.h (rev 0)
+++ grass-addons/grass7/raster/r.sim.water.mp/simlib.h 2017-04-28 04:43:41 UTC (rev 70973)
@@ -0,0 +1,123 @@
+#ifndef __SIMLIB_MP_H__
+#define __SIMLIB_MP_H__
+
+/*! \file simlib.h
+ * \brief This is the interface for the simlib (SIMWE) library.
+ */
+
+#define PARALLEL
+
+#ifdef PARALLEL
+#include <omp.h>
+#define NUM_THREADS "4"
+#endif
+
+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);
+
+#ifndef PARALLEL
+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;
+};
+#endif
+
+#ifdef PARALLEL
+
+int threads;
+
+time_t timer;
+
+void printTimeDiff(const char* message);
+
+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, *threads;
+};
+#endif
+
+struct flags
+{
+ struct Flag *mscale, *tserie;
+};
+
+#endif /* __SIMLIB_MP_H__ */
Added: grass-addons/grass7/raster/r.sim.water.mp/spearfish.sh
===================================================================
--- grass-addons/grass7/raster/r.sim.water.mp/spearfish.sh (rev 0)
+++ grass-addons/grass7/raster/r.sim.water.mp/spearfish.sh 2017-04-28 04:43:41 UTC (rev 70973)
@@ -0,0 +1,61 @@
+#!/bin/sh
+
+# sample script for Spearfish
+# written by
+# andreas.philipp geo.uni-augsburg.de
+# modified by JH, HM, MN, SG
+
+dem=elevation.10m
+output=simwe
+g.region raster=${output}
+g.region n=4920800 s=4917800 w=602500 e=606000
+g.region -p
+
+# rainfall [mm/hr]
+rain01=25
+
+# infitration [mm/hr]
+infil0=0.
+
+echo "Preparing input maps..."
+r.slope.aspect --o elevation=$dem dx=${output}_dx dy=${output}_dy
+r.mapcalc --o expr="${output}_rain =if(${dem}, $rain01, null())"
+r.mapcalc --o expr="${output}_manin=if(${dem}, soils * 0.005, null())"
+r.mapcalc --o expr="${output}_infil=if(${dem}, $infil0, null())"
+r.mapcalc --o expr="${output}_null=if(${dem}, float(0.0), null())"
+
+# Create the observation points
+v.random --o output=observation_points n=5 -z
+
+echo "r.sim.water --o -t elevation=${dem} dx=${output}_dx dy=${output}_dy rain=${output}_rain man=${output}_manin infil=${output}_infil depth=${output}_depth disch=${output}_disch err=${output}_err outwalk=${output}_walker niter=5 outiter=1 observation=observation_points logfile=simwe.log"
+r.sim.water --o -t elevation=${dem} dx=${output}_dx dy=${output}_dy \
+ rain=${output}_rain man=${output}_manin infil=${output}_infil \
+ depth=${output}_depth disch=${output}_disch err=${output}_err \
+ outwalk=${output}_walker niter=5 outiter=1 observation=observation_points \
+ logfile=simwe.log
+
+echo "Plotting the observation points with R. Result in Rplots.pdf"
+R --vanilla -e "df = read.table(\"simwe.log\", header=1); par(mfrow=c(3,2)); plot(df\$STEP, df\$CAT0001, type=\"l\"); plot(df\$STEP, df\$CAT0002, type=\"l\"); plot(df\$STEP, df\$CAT0003, type=\"l\"); plot(df\$STEP, df\$CAT0004, type=\"l\");plot(df\$STEP, df\$CAT0005, type=\"l\")"
+
+echo "Export of manning, soil and gradients"
+r.out.vtk elevation=${dem} input=${output}_manin,soils vectormaps=${output}_dx,${output}_dy,${output}_null output=manning_soils_gradient.vtk null=0.0
+
+echo "Build topology and exporting walker vector points for each time step"
+for i in `g.list vect | grep ${output}` ; do
+ v.build $i
+ echo "v.out.vtk input=$i output=$i.vtk"
+ v.out.vtk input=$i output=$i.vtk
+done
+
+echo "Exporting the raster maps for each time step"
+for i in `g.list rast | grep ${output}` ; do
+ echo "r.out.vtk elevation=${dem} input=$i output=$i.vtk null=0.0"
+ r.out.vtk elevation=${dem} input=$i output=$i.vtk null=0.0
+done
+
+echo "Now open paraview from this command line and load the vtk files as time series using the File->Open menu entry"
+echo "Step throu the time steps and adjust the color tables"
+
+# cleanup
+g.remove --q -f type=raster name=${output}_dx,${output}_dy,${output}_rain,${output}_manin,${output}_infil,${output}_null
+g.remove --q -f type=vector name=observation_points
Property changes on: grass-addons/grass7/raster/r.sim.water.mp/spearfish.sh
___________________________________________________________________
Added: svn:executable
+ *
Added: grass-addons/grass7/raster/r.sim.water.mp/utils.c
===================================================================
--- grass-addons/grass7/raster/r.sim.water.mp/utils.c (rev 0)
+++ grass-addons/grass7/raster/r.sim.water.mp/utils.c 2017-04-28 04:43:41 UTC (rev 70973)
@@ -0,0 +1,60 @@
+/* output.c (simlib), 20.nov.2002, JH */
+
+#include "waterglobs.h"
+
+double amax1(double arg1, double arg2)
+{
+ double res;
+
+ if (arg1 >= arg2) {
+ res = arg1;
+ }
+ else {
+ res = arg2;
+ }
+
+ return res;
+}
+
+double amin1(double arg1, double arg2)
+{
+ double res;
+
+ if (arg1 <= arg2) {
+ res = arg1;
+ }
+ else {
+ res = arg2;
+ }
+
+ return res;
+}
+
+int min(int arg1, int arg2)
+{
+ int res;
+
+ if (arg1 <= arg2) {
+ res = arg1;
+ }
+ else {
+ res = arg2;
+ }
+
+ return res;
+}
+
+int max(int arg1, int arg2)
+{
+ int res;
+
+ if (arg1 >= arg2) {
+ res = arg1;
+ }
+ else {
+ res = arg2;
+ }
+
+ return res;
+}
+
Added: grass-addons/grass7/raster/r.sim.water.mp/waterglobs.h
===================================================================
--- grass-addons/grass7/raster/r.sim.water.mp/waterglobs.h (rev 0)
+++ grass-addons/grass7/raster/r.sim.water.mp/waterglobs.h 2017-04-28 04:43:41 UTC (rev 70973)
@@ -0,0 +1,112 @@
+#ifndef __WATERGLOBS_MP_H__
+#define __WATERGLOBS_MP_H__
+
+#define EPS 1.e-7
+#define MAXW 20000000
+#define UNDEF -9999
+
+#include <grass/raster.h>
+
+extern char *elevin;
+extern char *dxin;
+extern char *dyin;
+extern char *rain;
+extern char *infil;
+extern char *traps;
+extern char *manin;
+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;
+
+extern char *wdepth;
+extern char *detin;
+extern char *tranin;
+extern char *tauin;
+extern char *tc;
+extern char *et;
+extern char *conc;
+extern char *flux;
+extern char *erdep;
+
+extern char *rainval;
+extern char *maninval;
+extern char *infilval;
+
+struct seed
+{
+ long int is1, is2;
+};
+
+extern struct seed seed;
+
+struct _points
+{
+ double *x; /* x coor for each point */
+ double *y; /* y coor for each point*/
+ int *cats; /* Category for each point */
+ int npoints; /* Number of observation points */
+ int npoints_alloc; /* Number of allocated points */
+ FILE *output; /* Output file descriptor */
+ int is_open; /* Set to 1 if open, 0 if closed */
+};
+
+extern struct _points points;
+extern void erod(double **);
+extern int output_et(void);
+extern double simwe_rand(void);
+extern double gasdev(void);
+extern void gasdev_for_paralel(double *x, double *y);
+extern double amax1(double, double);
+extern double amin1(double, double);
+extern int min(int, int);
+extern int max(int, int);
+extern void create_observation_points();
+
+extern double xmin, ymin, xmax, ymax;
+extern double mayy, miyy, maxx, mixx;
+extern int mx, my;
+extern int mx2, my2;
+
+extern double bxmi, bymi, bxma, byma, bresx, bresy;
+extern int maxwab;
+extern double step, conv;
+
+extern double frac;
+extern double bxmi, bymi;
+
+extern float **zz, **cchez;
+extern double **v1, **v2, **slope;
+extern double **gama, **gammas, **si, **inf, **sigma;
+extern float **dc, **tau, **er, **ct, **trap;
+extern float **dif;
+
+extern double vavg[MAXW][2], stack[MAXW][3], w[MAXW][3];
+extern int iflag[MAXW];
+
+extern double hbeta;
+extern double hhmax, sisum, vmean;
+extern double infsum, infmean;
+extern int maxw, maxwa, nwalk;
+extern double rwalk, bresx, bresy, xrand, yrand;
+extern double stepx, stepy, xp0, yp0;
+extern double chmean, si0, deltap, deldif, cch, hhc, halpha;
+extern double eps;
+extern int maxwab, nstack;
+extern int iterout, mx2o, my2o;
+extern int miter, nwalka;
+extern double timec;
+extern int ts, timesec;
+
+extern double rain_val;
+extern double manin_val;
+extern double infil_val;
+
+extern struct History history; /* holds meta-data (title, comments,..) */
+
+#endif /* __WATERGLOBS_MP_H__ */
More information about the grass-commit
mailing list