Skip to content

Commit

Permalink
r.sim: cleanup of global variables (#5044)
Browse files Browse the repository at this point in the history 
* group global variables into structures
* remove unused code
* limit scope of some variables
* fix small bug in parallelization code (`nwalka` variable)
  • Loading branch information
@petrasovaa
petrasovaa authored Feb 11, 2025
1 parent be43b74 commit a27e3c1
Show file tree
Hide file tree
Showing 8 changed files with 441 additions and 638 deletions.
110 changes: 46 additions & 64 deletions raster/r.sim/r.sim.sediment/main.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -321,50 +321,32 @@ int main(int argc, char *argv[])

G_get_set_window(&cellhd);

Geometry geometry = {0};
Settings settings = {0};
Setup setup = {0};
Simulation sim = {0};
settings.hhmax = settings.halpha = settings.hbeta = 0;
settings.ts = false;

WaterParams_init(&wp);

wp.conv = G_database_units_to_meters_factor();

wp.mixx = cellhd.west * wp.conv;
wp.maxx = cellhd.east * wp.conv;
wp.miyy = cellhd.south * wp.conv;
wp.mayy = cellhd.north * wp.conv;

wp.stepx = cellhd.ew_res * wp.conv;
wp.stepy = cellhd.ns_res * wp.conv;
/* wp.step = amin1(wp.stepx,wp.stepy); */
wp.step = (wp.stepx + wp.stepy) / 2.;
wp.mx = cellhd.cols;
wp.my = cellhd.rows;
wp.xmin = 0.;
wp.ymin = 0.;
wp.xp0 = wp.xmin + wp.stepx / 2.;
wp.yp0 = wp.ymin + wp.stepy / 2.;
wp.xmax = wp.xmin + wp.stepx * (float)wp.mx;
wp.ymax = wp.ymin + wp.stepy * (float)wp.my;
wp.hhc = wp.hhmax = 0.;

#if 0
wp.bxmi = 2093113. * wp.conv;
wp.bymi = 731331. * wp.conv;
wp.bxma = 2093461. * wp.conv;
wp.byma = 731529. * wp.conv;
wp.bresx = 2. * wp.conv;
wp.bresy = 2. * wp.conv;
wp.maxwab = 100000;

wp.mx2o = (int)((wp.bxma - wp.bxmi) / wp.bresx);
wp.my2o = (int)((wp.byma - wp.bymi) / wp.bresy);

/* relative small box coordinates: leave 1 grid layer for overlap */

wp.bxmi = wp.bxmi - wp.mixx + wp.stepx;
wp.bymi = wp.bymi - wp.miyy + wp.stepy;
wp.bxma = wp.bxma - wp.mixx - wp.stepx;
wp.byma = wp.byma - wp.miyy - wp.stepy;
wp.mx2 = wp.mx2o - 2 * ((int)(wp.stepx / wp.bresx));
wp.my2 = wp.my2o - 2 * ((int)(wp.stepy / wp.bresy));
#endif
geometry.conv = G_database_units_to_meters_factor();

geometry.mixx = cellhd.west * geometry.conv;
geometry.miyy = cellhd.south * geometry.conv;

geometry.stepx = cellhd.ew_res * geometry.conv;
geometry.stepy = cellhd.ns_res * geometry.conv;
/* geometry.step = amin1(geometry.stepx,geometry.stepy); */
geometry.step = (geometry.stepx + geometry.stepy) / 2.;
geometry.mx = cellhd.cols;
geometry.my = cellhd.rows;
geometry.xmin = 0.;
geometry.ymin = 0.;
geometry.xp0 = geometry.xmin + geometry.stepx / 2.;
geometry.yp0 = geometry.ymin + geometry.stepy / 2.;
geometry.xmax = geometry.xmin + geometry.stepx * (float)geometry.mx;
geometry.ymax = geometry.ymin + geometry.stepy * (float)geometry.my;

wp.elevin = parm.elevin->answer;
wp.wdepth = parm.wdepth->answer;
Expand Down Expand Up @@ -400,35 +382,35 @@ int main(int argc, char *argv[])
G_message(_("Number of threads: %d"), threads);

/* sscanf(parm.nwalk->answer, "%d", &wp.maxwa); */
sscanf(parm.niter->answer, "%d", &wp.timesec);
sscanf(parm.mintimestep->answer, "%lf", &wp.mintimestep);
sscanf(parm.outiter->answer, "%d", &wp.iterout);
sscanf(parm.niter->answer, "%d", &settings.timesec);
sscanf(parm.outiter->answer, "%d", &settings.iterout);
sscanf(parm.mintimestep->answer, "%lf", &settings.mintimestep);
/* sscanf(parm.density->answer, "%d", &wp.ldemo); */
sscanf(parm.diffc->answer, "%lf", &wp.frac);
sscanf(parm.diffc->answer, "%lf", &settings.frac);
sscanf(parm.maninval->answer, "%lf", &wp.manin_val);

/* Recompute timesec from user input in minutes
* to real timesec in seconds */
wp.timesec = wp.timesec * 60;
wp.iterout = wp.iterout * 60;
if ((wp.timesec / wp.iterout) > 100)
settings.timesec = settings.timesec * 60;
settings.iterout = settings.iterout * 60;
if ((settings.timesec / settings.iterout) > 100)
G_message(_("More than 100 files are going to be created !!!!!"));

/* compute how big the raster is and set this to appr 2 walkers per cell */
if (parm.nwalk->answer == NULL) {
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);
sim.maxwa = geometry.mx * geometry.my * 2;
sim.rwalk = (double)(geometry.mx * geometry.my * 2.);
G_message(_("default nwalk=%d, rwalk=%f"), sim.maxwa, sim.rwalk);
}
else {
sscanf(parm.nwalk->answer, "%d", &wp.maxwa);
wp.rwalk = (double)wp.maxwa;
sscanf(parm.nwalk->answer, "%d", &sim.maxwa);
sim.rwalk = (double)sim.maxwa;
}
/*rwalk = (double) maxwa; */

if (wp.conv != 1.0)
G_message(_("Using metric conversion factor %f, step=%f"), wp.conv,
wp.step);
if (geometry.conv != 1.0)
G_message(_("Using metric conversion factor %f, step=%f"),
geometry.conv, geometry.step);

wp.observation = parm.observation->answer;
wp.logfile = parm.logfile->answer;
Expand All @@ -437,24 +419,24 @@ int main(int argc, char *argv[])
if ((wp.tc == NULL) && (wp.et == NULL) && (wp.conc == NULL) &&
(wp.flux == NULL) && (wp.erdep == NULL))
G_warning(_("You are not outputting any raster or site files"));
ret_val = input_data();
ret_val = input_data(geometry.my, geometry.mx, &sim);
if (ret_val != 1)
G_fatal_error(_("Input failed"));

alloc_grids_sediment();
alloc_grids_sediment(&geometry);

grad_check();
init_grids_sediment();
grad_check(&setup, &geometry, &settings);
init_grids_sediment(&setup, &geometry);
/* treba dat output pre topoerdep */
main_loop();
main_loop(&setup, &geometry, &settings, &sim);

/* always true for sediment? */
if (wp.tserie == NULL) {
ii = output_data(0, 1.);
ii = output_data(0, 1., &setup, &geometry, &settings, &sim);
if (ii != 1)
G_fatal_error(_("Cannot write raster maps"));
}
free_walkers();
free_walkers(&sim);

/* Exit with Success */
exit(EXIT_SUCCESS);
Expand Down
87 changes: 45 additions & 42 deletions raster/r.sim/r.sim.water/main.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -355,35 +355,38 @@ int main(int argc, char *argv[])

G_get_set_window(&cellhd);

Geometry geometry = {0};
Settings settings = {0};
Setup setup = {0};
Simulation sim = {0};

WaterParams_init(&wp);

wp.conv = G_database_units_to_meters_factor();
geometry.conv = G_database_units_to_meters_factor();

G_debug(3, "Conversion factor is set to: %f", wp.conv);
G_debug(3, "Conversion factor is set to: %f", geometry.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;
geometry.mixx = geometry.conv * cellhd.west;
geometry.miyy = geometry.conv * cellhd.south;

wp.stepx = cellhd.ew_res * wp.conv;
wp.stepy = cellhd.ns_res * wp.conv;
geometry.stepx = cellhd.ew_res * geometry.conv;
geometry.stepy = cellhd.ns_res * geometry.conv;
/* step = amin1(stepx,stepy); */
wp.step = (wp.stepx + wp.stepy) / 2.;
wp.mx = cellhd.cols;
wp.my = cellhd.rows;
geometry.step = (geometry.stepx + geometry.stepy) / 2.;
geometry.mx = cellhd.cols;
geometry.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;
geometry.xmin = 0.;
geometry.ymin = 0.;
geometry.xp0 = geometry.xmin + geometry.stepx / 2.;
geometry.yp0 = geometry.ymin + geometry.stepy / 2.;
geometry.xmax = geometry.xmin + geometry.stepx * (float)geometry.mx;
geometry.ymax = geometry.ymin + geometry.stepy * (float)geometry.my;

G_debug(3, "xmax: %f, ymax: %f", wp.xmax, wp.ymax);
G_debug(3, "xmax: %f, ymax: %f", geometry.xmax, geometry.ymax);

wp.ts = flag.tserie->answer;
settings.ts = flag.tserie->answer;

wp.elevin = parm.elevin->answer;
wp.dxin = parm.dxin->answer;
Expand All @@ -399,13 +402,13 @@ int main(int argc, char *argv[])

G_debug(3, "Parsing numeric parameters");

sscanf(parm.niter->answer, "%d", &wp.timesec);
sscanf(parm.mintimestep->answer, "%lf", &wp.mintimestep);
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);
sscanf(parm.niter->answer, "%d", &settings.timesec);
sscanf(parm.outiter->answer, "%d", &settings.iterout);
sscanf(parm.mintimestep->answer, "%lf", &settings.mintimestep);
sscanf(parm.diffc->answer, "%lf", &settings.frac);
sscanf(parm.hmax->answer, "%lf", &settings.hhmax);
sscanf(parm.halpha->answer, "%lf", &settings.halpha);
sscanf(parm.hbeta->answer, "%lf", &settings.hbeta);

G_debug(3, "Parsing rain parameters");

Expand Down Expand Up @@ -516,43 +519,43 @@ int main(int argc, char *argv[])

/* 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)
settings.timesec = settings.timesec * 60.0;
settings.iterout = settings.iterout * 60.0;
if ((settings.timesec / settings.iterout) > 100.0 && settings.ts)
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);
sim.maxwa = geometry.mx * geometry.my * 2;
sim.rwalk = (double)(geometry.mx * geometry.my * 2.);
G_message(_("default nwalk=%d, rwalk=%f"), sim.maxwa, sim.rwalk);
}
else {
sscanf(parm.nwalk->answer, "%d", &wp.maxwa);
wp.rwalk = (double)wp.maxwa;
sscanf(parm.nwalk->answer, "%d", &sim.maxwa);
sim.rwalk = (double)sim.maxwa;
}

/* rwalk = (double) maxwa; */

if (wp.conv != 1.0)
G_message(_("Using metric conversion factor %f, step=%f"), wp.conv,
wp.step);
if (geometry.conv != 1.0)
G_message(_("Using metric conversion factor %f, step=%f"),
geometry.conv, geometry.step);

wp.observation = parm.observation->answer;
wp.logfile = parm.logfile->answer;
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();
ret_val = input_data(geometry.my, geometry.mx, &sim);
if (ret_val != 1)
G_fatal_error(_("Input failed"));

alloc_grids_water();
alloc_grids_water(&geometry);

grad_check();
main_loop();
free_walkers();
grad_check(&setup, &geometry, &settings);
main_loop(&setup, &geometry, &settings, &sim);
free_walkers(&sim);

/* Exit with Success */
exit(EXIT_SUCCESS);
Expand Down
25 changes: 15 additions & 10 deletions raster/r.sim/simlib/erod.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,41 +22,46 @@
#include <grass/linkm.h>

#include <grass/waterglobs.h>
#include <grass/simlib.h>

/* divergence computation from a given field */

void erod(double **hw)
void erod(double **hw, const Setup *setup, const Geometry *geometry)
{
/* 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++) {
for (k = 0; k < geometry->my; k++) {
for (l = 0; l < geometry->mx; l++) {
lp = max(0, l - 2);
l1 = lp + 1;
kp = max(0, k - 2);
k1 = kp + 1;
ln = min(mx - 1, l + 1);
ln = min(geometry->mx - 1, l + 1);
l2 = ln - 1;
kn = min(my - 1, k + 1);
kn = min(geometry->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;
dxp = (v1[k][lp] * hw[k][lp] - v1[k][l1] * hw[k][l1]) /
geometry->stepx;
dxn = (v1[k][l2] * hw[k][l2] - v1[k][ln] * hw[k][ln]) /
geometry->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;
dyp = (v2[kp][l] * hw[kp][l] - v2[k1][l] * hw[k1][l]) /
geometry->stepy;
dyn = (v2[k2][l] * hw[k2][l] - v2[kn][l] * hw[kn][l]) /
geometry->stepy;
dya = 0.5 * (dyp + dyn);

er[k][l] = (dxa + dya) / deltap;
er[k][l] = (dxa + dya) / setup->deltap;
}
else
er[k][l] = UNDEF;
Expand Down
Loading

0 comments on commit a27e3c1

Please sign in to comment.