r.sim: add parallelization support to r.sim.water and r.sim.sediment - ported r.sim.water.mp from addons by Jaroslav Hofierka et al

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

raster/r.sim/r.sim.sediment/Makefile

@@ -4,10 +4,10 @@ PGM=r.sim.sediment
 
 EXTRA_CLEAN_DIRS=doxygenhtml
 
-LIBES     = $(SIMLIB) $(GMATHLIB) $(GISLIB)
+LIBES     = $(SIMLIB) $(GMATHLIB) $(GISLIB) $(OMPLIB)
 DEPENDENCIES = $(SIMDEP) $(GMATHDEP) $(GISDEP)
 EXTRA_INC = $(VECT_INC)
-EXTRA_CFLAGS = -I ../simlib $(VECT_CFLAGS) 
+EXTRA_CFLAGS = -I ../simlib $(VECT_CFLAGS) $(OMPCFLAGS)
 
 include $(MODULE_TOPDIR)/include/Make/Module.make
 

raster/r.sim/r.sim.sediment/main.c

@@ -100,6 +100,7 @@ char msg[1024];
 int main(int argc, char *argv[])
 {
     int ii;
+    int threads;
     int ret_val;
     struct Cell_head cellhd;
     struct WaterParams wp;
@@ -278,6 +279,15 @@ int main(int argc, char *argv[])
         _("Automatically generates random seed for random number"
           " generator (use when you don't want to provide the seed option)");
 
+    parm.threads = G_define_option();
+    parm.threads->key = "nprocs";
+    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");
+
     if (G_parser(argc, argv))
 	exit(EXIT_FAILURE);
 
@@ -358,6 +368,16 @@ int main(int argc, char *argv[])
     wp.erdep = parm.erdep->answer;
     wp.outwalk = parm.outwalk->answer; 
 
+    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 threads: %d"), threads);
+
     /*      sscanf(parm.nwalk->answer, "%d", &wp.maxwa); */
     sscanf(parm.niter->answer, "%d", &wp.timesec);
     sscanf(parm.outiter->answer, "%d", &wp.iterout);

raster/r.sim/r.sim.water/Makefile

@@ -4,10 +4,10 @@ PGM=r.sim.water
 
 EXTRA_CLEAN_DIRS=doxygenhtml
 
-LIBES     = $(SIMLIB) $(GMATHLIB) $(GISLIB)
+LIBES     = $(SIMLIB) $(GMATHLIB) $(GISLIB) $(OMPLIB)
 DEPENDENCIES = $(SIMDEP) $(GMATHDEP) $(GISDEP)
 EXTRA_INC = $(VECT_INC)
-EXTRA_CFLAGS = -I ../simlib $(VECT_CFLAGS) 
+EXTRA_CFLAGS = -I ../simlib $(VECT_CFLAGS) $(OMPCFLAGS)
 
 include $(MODULE_TOPDIR)/include/Make/Module.make
 

raster/r.sim/r.sim.water/main.c

@@ -90,7 +90,6 @@
 
 #include <grass/simlib.h>
 
-
 /****************************************/
 /* MAIN                                 */
 
@@ -98,6 +97,7 @@
 int main(int argc, char *argv[])
 {
     int ii;
+    int threads;
     int ret_val;
     double x_orig, y_orig;
     struct GModule *module;
@@ -314,6 +314,15 @@ int main(int argc, char *argv[])
         _("Automatically generates random seed for random number"
           " generator (use when you don't want to provide the seed option)");
 
+     parm.threads = G_define_option();
+     parm.threads->key = "nprocs";
+     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");
+
     if (G_parser(argc, argv))
 	exit(EXIT_FAILURE);
 
@@ -387,6 +396,16 @@ int main(int argc, char *argv[])
 
     G_debug(3, "Parsing rain parameters");
 
+    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 threads: %d"), threads);
+
     /* if no rain map input, then: */
     if (parm.rain->answer == NULL) {
 	/*Check for Rain Unique Value Input */

raster/r.sim/simlib/Makefile

@@ -1,9 +1,9 @@
 MODULE_TOPDIR = ../../..
 
-LIBES = $(GISLIB) $(DATETIMELIB) $(GMATHLIB) $(RASTERLIB) $(VECTORLIB)
+LIBES = $(GISLIB) $(DATETIMELIB) $(GMATHLIB) $(RASTERLIB) $(VECTORLIB) $(OMPLIB)
 DEPENDENCIES = $(GISDEP) $(BITMAPDEP) $(DBMIDEP) $(GMATHDEP) $(LINKMDEP) $(XDRDEP) $(SITESDEP) $(VECTORDEP) $(RASTERDEP) $(VECTORDEP)
 EXTRA_INC = $(VECT_INC)
-EXTRA_CFLAGS = $(VECT_CFLAGS)
+EXTRA_CFLAGS = $(VECT_CFLAGS) $(OMPCFLAGS)
 
 LIB = SIM
 

raster/r.sim/simlib/hydro.c

@@ -268,7 +268,14 @@ void main_loop(void)
 	    nwalka = 0;
 	    nstack = 0;
 
-	    for (lw = 0; lw < nwalk; lw++) {
+#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;
+
+	    for (lw = min_loop; lw < max_loop; lw++) {
 		if (w[lw][2] > EPS) {	/* check the walker weight */
 		    ++nwalka;
 		    l = (int)((w[lw][0] + stxm) / stepx) - mx - 1;
@@ -305,6 +312,7 @@ void main_loop(void)
 			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 */
@@ -329,9 +337,6 @@ void main_loop(void)
 			    }
 			}
 
-			gaux = gasdev();
-			gauy = gasdev();
-
 			w[lw][0] += (velx + dif[k][l] * gaux);	/* move the walker */
 			w[lw][1] += (vely + dif[k][l] * gauy);
 
@@ -358,7 +363,7 @@ void main_loop(void)
 		    }
 		}
             } /* 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 */

raster/r.sim/simlib/random.c

@@ -45,3 +45,17 @@ double gasdev(void)
     }
     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;
+}

raster/r.sim/simlib/simlib.h

@@ -5,6 +5,12 @@
  * \brief This is the interface for the simlib (SIMWE) library.
  */
 
+
+#define NUM_THREADS "1"
+#if defined(_OPENMP)
+#include <omp.h>
+#endif
+
 struct WaterParams
 {
     double xmin, ymin, xmax, ymax;
@@ -85,7 +91,7 @@ struct options
 	*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, *seed;
+	*infilval, *logfile, *seed, *threads;
 };
 
 struct flags

raster/r.sim/simlib/waterglobs.h

@@ -61,6 +61,7 @@ extern void erod(double **);
 extern int output_et(void);
 extern double simwe_rand(void);
 extern double gasdev(void);
+extern void gasdev_for_paralel(double *, double *);
 extern double amax1(double, double);
 extern double amin1(double, double);
 extern int min(int, int);