From 9ac2a068af1ab2ed906c7af7ac438bf35143ad1e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Edouard=20Choini=C3=A8re?= <27212526+echoix@users.noreply.github.com> Date: Mon, 3 Jan 2022 14:36:29 -0500 Subject: [PATCH 01/17] Add english comments --- src/raster/r.damflood/main.c | 21 ++++++++++++--------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/src/raster/r.damflood/main.c b/src/raster/r.damflood/main.c index da61ad9bb7..812eeeea85 100644 --- a/src/raster/r.damflood/main.c +++ b/src/raster/r.damflood/main.c @@ -53,6 +53,7 @@ //#define hmin 0.01 /* lo calcolo dopo in funzione della velocita' massima + (Calculated later as a function of maximum velocity) #define timestep 0.1 */ /* @@ -731,7 +732,7 @@ int main(int argc, char *argv[]){ if (method==1 || method==2) { num_cell=0; - /* cerco il lago */ + /* Search for lakes (cerco il lago) */ for (row = 0; row < nrows; row++) { for (col = 0; col < ncols; col++) @@ -750,7 +751,7 @@ int main(int argc, char *argv[]){ { for (col = 0; col < ncols; col++) { - /* rottura diga */ + /* Dam break (rottura diga) */ if (m_DAMBREAK[row][col] > 0){ num_break++; G_message("(%d,%d)Cell Dam Breach n° %d",row,col,num_break); @@ -781,6 +782,7 @@ int main(int argc, char *argv[]){ /**************************************/ /* timestep in funzione di V_0 e res */ + /* timestep as a function of V_0 and res */ /**************************************/ //timestep=0.01; //timestep= ((res_ns+res_ew)/2.0)/(vel_max*50.0); @@ -793,7 +795,7 @@ int main(int argc, char *argv[]){ // Uniform drop in of lake (method=1 or method =2) //***************************************************************************** if (method==1 || method==2){ - /* calcolo l'abbassamento sul lago*/ + /* calcolo l'abbassamento sul lago (Calculation of the lowering of the lake) */ if (num_cell!=0) { fall = (volume) / (num_cell * res_ew * res_ns); //printf("volume=%f, fall=%f\n",volume,fall); @@ -805,12 +807,13 @@ int main(int argc, char *argv[]){ for (col = 1; col < ncols-1; col++) { if (m_DAMBREAK[row][col]>0){ - // ragiona se ha senso + // ragiona se ha senso (I think it makes sense) m_h2[row][col]=m_h1[row][col]-fall; if (m_h2[row][col]<=0) { m_h2[row][col]=0.0; if (m_h1[row][col]>0) { // questo warning va modificato perchè vale per ogni cella ---> bisogna metterne uno generico che valga quando tutte le celle sono con h=0 + // (This warning must be modified since it is valid for every cell) ---> (You have to put a generic one that is valid when all cells have h=0) num_break--; if (num_break==0){ if (warn1==0){ @@ -837,7 +840,7 @@ int main(int argc, char *argv[]){ { for (col = 0; col < ncols; col++) { - /* rottura diga */ + /* Dam break (rottura diga) */ if (m_DAMBREAK[row][col] > 0){ m_z[row][col]=m_z[row][col]-m_DAMBREAK[row][col]; m_DAMBREAK[row][col]=-1.0; @@ -961,10 +964,10 @@ int main(int argc, char *argv[]){ /* qualcosa non va in questo if */ /*----------------------------------------------------------------------------------------------*/ - /* controllo se devo scrivere outputs */ + /* controllo se devo scrivere outputs (Check if we need to write outputs) */ if ((input_DELTAT->answer != NULL)||(parm.opt_t->answer != NULL)) { if ((((m*DELTAT-t) <= timestep && (m*DELTAT) < TSTOP) && (input_DELTAT->answer != NULL)) || ((pp < ntimes && (times[pp]-t) < timestep) && (parm.opt_t->answer != NULL))) { - if ((m*DELTAT-t) <= timestep && m*DELTAT < TSTOP) { /* devo cambiare il nome del raster e aggiungere ogni volta _timestep*/ + if ((m*DELTAT-t) <= timestep && m*DELTAT < TSTOP) { /* devo cambiare il nome del raster e aggiungere ogni volta _timestep (We need to change the raster name and add _timestep at each time)*/ if (OUT_H) { sprintf(name1,"%s%d",OUT_H,m*DELTAT); } @@ -1041,7 +1044,7 @@ int main(int argc, char *argv[]){ } } /* end_col*/ - /*copia righe !!! */ + /*copia righe (Copy lines)!!! */ if (OUT_VEL) { Rast_put_d_row(outfd_VEL,outrast_VEL); } @@ -1318,7 +1321,7 @@ for (row = 0; row < nrows; row++){ } //G_message("pippo, row=%d e nrows=%d", row, nrows); } // end row -/* chiudi i file */ +/* chiudi i file (Close files) */ if (OUT_H) { From 7eade51e291adcbf977fa1a73ff54a6fc5134393 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Edouard=20Choini=C3=A8re?= <27212526+echoix@users.noreply.github.com> Date: Mon, 3 Jan 2022 14:35:04 -0500 Subject: [PATCH 02/17] Spelling fixes in r.damflood --- src/raster/r.damflood/main.c | 50 ++++++++++++++++++------------------ 1 file changed, 25 insertions(+), 25 deletions(-) diff --git a/src/raster/r.damflood/main.c b/src/raster/r.damflood/main.c index 812eeeea85..e94bc64be7 100644 --- a/src/raster/r.damflood/main.c +++ b/src/raster/r.damflood/main.c @@ -33,7 +33,7 @@ #include #include /* function here defined */ -#include "SWE.h" /*function that solve the shallow water equations*/ +#include "SWE.h" /*function that solves the shallow water equations*/ //#include @@ -166,7 +166,7 @@ int main(int argc, char *argv[]){ startpt *a_start,*tmp_start;*/ struct Cell_head cellhd; /* it stores region information and header information of rasters */ struct Cell_head window; - //struct History history; /* holds meta-data */ + //struct History history; /* holds metadata */ /* input-output raster file descriptor */ int infd_ELEV,infd_LAKE,infd_DAMBREAK, infd_MANNING, infd_U, infd_V; @@ -544,7 +544,7 @@ int main(int argc, char *argv[]){ G_fatal_error("You choose flow direction map without flow velocity prefix name"); } - /* check legal output name */ + /* Check legal output name */ if (OUT_H) { if (G_legal_filename (OUT_H) < 0) G_fatal_error (_("[%s] is an illegal name"), OUT_H); @@ -569,7 +569,7 @@ int main(int argc, char *argv[]){ if (G_legal_filename (OUT_WAVEFRONT) < 0) G_fatal_error (_("[%s] is an illegal name"), OUT_WAVEFRONT); } - /* type of dam failure*/ + /* Type of dam failure*/ if (strcmp(opt.met->answer, "uniform drop in of lake") == 0) method=1; else if (strcmp(opt.met->answer, "small dam breach") == 0) @@ -596,14 +596,14 @@ int main(int argc, char *argv[]){ } else if (input_V->answer != NULL && input_U->answer == NULL ) { inrast_V = Rast_allocate_d_buf(); } - /* get windows rows & cols */ + /* Get windows rows & cols */ nrows = Rast_window_rows(); ncols = Rast_window_cols(); G_get_window(&window); res_ew = window.ew_res; res_ns = window.ns_res; - /* allocate memory matrix */ + /* Allocate memory matrix */ m_DAMBREAK = G_alloc_fmatrix(nrows,ncols); m_m = G_alloc_fmatrix(nrows,ncols); m_z = G_alloc_fmatrix(nrows,ncols); @@ -639,7 +639,7 @@ int main(int argc, char *argv[]){ { G_percent (row, nrows, 2); - /* read a line input maps into buffers*/ + /* Read a line input maps into buffers*/ Rast_get_f_row (infd_ELEV, inrast_ELEV, row); Rast_get_d_row (infd_LAKE, inrast_LAKE, row); Rast_get_f_row (infd_DAMBREAK, inrast_DAMBREAK, row); @@ -661,10 +661,10 @@ int main(int argc, char *argv[]){ if (G_get_f_raster_row (infd_MANNING, inrast_MANNING, row) < 0) G_fatal_error (_("Could not read from <%s>"),MANNING);*/ - /* read every cell in the line buffers */ + /* Read every cell in the line buffers */ for (col = 0; col < ncols; col++) { - /* store values in memory matrix (attenzione valori nulli!)*/ + /* Store values in memory matrix (attenzione valori nulli!) (attention null values!)*/ m_DAMBREAK[row][col] = ((FCELL *) inrast_DAMBREAK)[col]; m_m[row][col] = ((FCELL *) inrast_MANNING)[col]; m_z[row][col] = ((FCELL *) inrast_ELEV)[col]; @@ -765,7 +765,7 @@ int main(int argc, char *argv[]){ if (vel_0 > vel_max) vel_max=vel_0; } else { - G_fatal_error(_("Don't find the dambreak - Please select a correct map or adjust the computational region")); + G_fatal_error(_("Didn't find the dambreak. Please select a correct map or adjust the computational region.")); } if (m_DAMBREAK[row][col] > 0) { @@ -778,7 +778,7 @@ int main(int argc, char *argv[]){ } } } - G_message("the number of lake cell is': %d\n", num_cell); + G_message("The number of lake cell is': %d\n", num_cell); /**************************************/ /* timestep in funzione di V_0 e res */ @@ -834,7 +834,7 @@ int main(int argc, char *argv[]){ } }}//end two for cicles } //end if - // there isn't interest to find where is the lake --> everywhere m_lake[row][col]=0 + // There isn't interest to find where is the lake --> everywhere m_lake[row][col]=0 if (method==3){ for (row = 0; row < nrows; row++) { @@ -854,7 +854,7 @@ int main(int argc, char *argv[]){ G_message("Model running"); - /* calculate time step loop */ + /* Calculate time step loop */ for(t=0; t<=TSTOP; t+=timestep){ //printf("************************************************\n"); //while(!getchar()){ } @@ -885,7 +885,7 @@ int main(int argc, char *argv[]){ if (reg_lim==0) { G_warning("At the time %.3f the computational region is smaller than inundation",t); reg_lim=1; - } /* warning message only a time */ + } /* Warning message only a time */ } /* velocities at the limit of computational region */ //******************************************************************** @@ -984,7 +984,7 @@ int main(int argc, char *argv[]){ sprintf(name3,"%s%s%s%d","opt_","dir_",OUT_VEL,pp); G_message("Time: %lf, writing an optional output maps %d",t, pp); } - /* controlling if we can write the raster */ + /* Controlling if we can write the raster */ if (OUT_H) { if ( (outfd_H = Rast_open_new (name1,DCELL_TYPE)) < 0) { G_fatal_error (_("Could not open <%s>"),name1); @@ -1000,7 +1000,7 @@ int main(int argc, char *argv[]){ G_fatal_error (_("Could not open <%s>"),name3); } } - /* allocate output buffer */ + /* Allocate output buffer */ if (OUT_VEL) { outrast_VEL = Rast_allocate_d_buf(); } @@ -1014,7 +1014,7 @@ int main(int argc, char *argv[]){ { for (col = 0; col < ncols; col++) { - /* copy matrix in buffer */ + /* Copy matrix in buffer */ if (OUT_VEL) { ((DCELL *) outrast_VEL)[col] = sqrt(pow(m_u1[row][col],2.0) + pow(m_v1[row][col],2.0)); } @@ -1055,7 +1055,7 @@ int main(int argc, char *argv[]){ Rast_put_d_row(outfd_VEL_DIR,outrast_VEL_DIR); } } /* end row */ - /* memory cleanup */ + /* Memory cleanup */ if (OUT_VEL) { G_free(outrast_VEL); } @@ -1065,7 +1065,7 @@ int main(int argc, char *argv[]){ if (flag_d->answer) { G_free(outrast_VEL_DIR); } - /* close the raster maps */ + /* Close the raster maps */ if (OUT_VEL) { Rast_close (outfd_VEL); } @@ -1090,7 +1090,7 @@ int main(int argc, char *argv[]){ //******************************************************************* -// write final flooding map +// Write final flooding map if(OUT_H) { sprintf(name1,"%s%d",OUT_H,TSTOP); } @@ -1151,7 +1151,7 @@ if (OUT_WAVEFRONT){ G_fatal_error (_("Could not open <%s>"),OUT_WAVEFRONT); } -/* allocate output buffer */ +/* Allocate output buffer */ if (OUT_H) { outrast_H = Rast_allocate_d_buf(); } @@ -1189,7 +1189,7 @@ G_percent(nrows, nrows, 1); /* finish it */ for (row = 0; row < nrows; row++){ G_percent (row, nrows, 2); for (col = 0; col < ncols; col++) { - /* copy matrix in buffer */ + /* Copy matrix in buffer */ if (OUT_H) { ((DCELL *) outrast_H)[col] = m_h1[row][col]; } @@ -1343,7 +1343,7 @@ for (row = 0; row < nrows; row++){ G_message("Writing the output final map %s", OUT_WAVEFRONT); } -G_percent(nrows, nrows, 1); /* finish it */ +G_percent(nrows, nrows, 1); /* Finish it */ if (OUT_VEL) { G_free(outrast_VEL); Rast_close (outfd_VEL); @@ -1389,9 +1389,9 @@ if (OUT_WAVEFRONT){ } //************************************************************************ -// da sistemare +// da sistemare (TODO/To fix) //************************************************************************ -/* add command line incantation to history file */ +/* Add command line incantation to history file */ //G_short_history(result, "raster", &history); //G_command_history(&history); //G_write_history(result, &history); From d80bfd3c85337b9a426244abc6a283f4a12d390d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Edouard=20Choini=C3=A8re?= <27212526+echoix@users.noreply.github.com> Date: Mon, 3 Jan 2022 18:24:35 -0500 Subject: [PATCH 03/17] r.damflood manual page revision --- src/raster/r.damflood/r.damflood.html | 241 ++++++++++++++++---------- 1 file changed, 149 insertions(+), 92 deletions(-) diff --git a/src/raster/r.damflood/r.damflood.html b/src/raster/r.damflood/r.damflood.html index b1a80450bf..78c54cc23d 100644 --- a/src/raster/r.damflood/r.damflood.html +++ b/src/raster/r.damflood/r.damflood.html @@ -1,124 +1,156 @@

DESCRIPTION

-r.damflood - The definition of flooding areas is of considerable importance for both the risk analysis and the emergency management. -This command, in particular, is an embedded GRASS GIS hydrodynamic 2D model that allows to obtain flooding area due to a failure -of a dam, given the geometry of the reservoir and of the downstream area, the initial conditions and the dam breach geometry. +r.damflood - The definition of flooding areas is of considerable +importance for both the risk analysis and the emergency management. +This module, in particular, is an embedded GRASS GIS hydrodynamic 2D model +that allows to obtain flooding areas due to a failure of a dam, given +the geometry of the reservoir and of the downstream area, +the initial conditions and the dam breach geometry.
-The numerical model solves the conservative form of the shallow water equations (SWE) using a finite volume method (FVM); -the intercell flux is computed by the "upwind method and the water-level gradient is evaluated by weighted -average of both upwind and downwind gradient. Additional details of the specific numerical scheme adopted in the model are +The numerical model solves the conservative form of the shallow water equations +(SWE) using a finite volume method (FVM); +the intercell flux is computed by the "upwind method and the water-level +gradient is evaluated by weighted average of both upwind and downwind gradient. +Additional details of the specific numerical scheme adopted in the model are presented in references [1].
-The command allows to generate raster time series, of water depth and flow velocity, with time resolution defined by user. -Each time series is identified by a number of raster maps named with a common prefix as specified by the user and the time instant -which it refers expressed in seconds from the dam failure, joined by the underscore character (e.g.; myvel_125, myvel_250, myvel_375, etc.).
- -Because this new module has been implemented with the aim to provide an instrument for risk assessment fully within a GIS environment, -it should be able to provide intensity maps directly applicable in those analyses.In floods, intensity generally corresponds -to the maximum flow depth, but in the particular case of flash floods, where velocities are normally high, -it is recommended to use as intensity indicator the maximum between the water depth and the product of water velocity and water depth. -For this reason, with this module, in addition to the water depth and velocity maps, the user can choose -a variety of output raster maps: maximum water depth, maximum water velocity, and maximum intensity raster maps.
- -In case on high numerical stability problem, the user is warned, and the simulation is stopped.
-

+The module allows to generate raster time series, of water depth and +flow velocity, with time resolution defined by the user. +Each time series is identified by a number of raster maps named with a common +prefix as specified by the user and the time instant which it refers +expressed in seconds from the dam failure, joined by the underscore character +(e.g.; myvel_125, myvel_250, myvel_375, etc.).
+ +Because this new module has been implemented with the aim to provide an +instrument for risk assessment fully within a GIS environment, it should be +able to provide intensity maps directly applicable in those analyses. +In floods, intensity generally corresponds to the maximum flow depth, +but in the particular case of flash floods, where velocities are normally high, +it is recommended to use as an intensity indicator the maximum between +the water depth and the product of water velocity and water depth. +For this reason, with this module, in addition to the water depth and velocity +maps, the user can choose a variety of output raster maps: maximum water depth, +maximum water velocity, and maximum intensity raster maps.
+ +In case of numerical stability problems, the user is warned, and the +simulation is stopped.
+
-Use -
-

-Requested input:
-The required input are:
-- a DTM including the lake bathimetry and the dam elevation over the ground [elev],
-- a map with the initial condition easily obtained with r.lake command [lake],
-- a dam breach width raster map [dambreak] which can be obtained using r.dam grass add-on script,
-- a Manning's roughness coefficient raster map, easily obtained from a reclassification of a land use map -(r.reclass) [manning],
-- the simulation time length expressed in seconds [tstop].

- - - - - -Output map and additional output options:
-First the user can set a specific time lag [deltat] expressed in seconds, that is used for the output map (depth and velocity) generation. -and also an additional series of instants [opt_t],expressed in seconds from the beginning of the simulation), used to generate further water flow depth and velocity maps -at desired precise times.
- -The user can choose between one of the following time series raster maps as output: - - flow depth [h],
- - flow velocity [vel],
- - - - a raster map with maximum water depth [hmax], relative flooding intensity [i_hmax], that is the product of water depth and velocity, and the relative time of occurence[t_hmax],
- - a raster map with maximum water velocity [vmax], relative flooding intensity [i_vmax], and the relative time of occurence[t_vmax],
- - a raster map with maximum flooding intensity [imax] and the relative time of occurence[t_imax].
- - a raster map with the time of arriving of the Wave-Front [wavefront]

-where and the raster maps are coded as "prefix" + "_" + "elapsed seconds": e.g. mydepth_125.

+

USAGE

+

Required inputs

+The required inputs are: +
    +
  • a DTM including the lake bathymetry and the dam elevation over + the ground [elev],
    • +
  • a map with the initial condition easily obtained with + r.lake + command [lake],
    • +
  • a dam breach width raster map [dambreak] which can be + obtained using r.dam grass add-on script,
    • +
  • a Manning's roughness coefficient raster map, easily obtained from a + reclassification of a land use map + (r.reclass) + [manning],
    • +
  • the simulation time length expressed in seconds + [tstop].
    • +
+ + +

Output map and additional output options

+First the user can set a specific time lag [deltat] expressed in +seconds, that is used for the output map (depth and velocity) +generation. +It is also possible to generate the output maps, like water flow depth and +velocity maps, at some precise times. These additional series of instants +are set with the opt_t parameter, expressed in seconds +from the beginning of the simulation. + +The user can choose between one of the following time series raster maps as +output: +
    +
  • flow depth [h],
    • +
  • flow velocity [vel],
    • + +
  • a raster map with maximum water depth [hmax], + relative flooding intensity [i_hmax], + that is the product of water depth and velocity, + and the relative time of occurrence [t_hmax],
    • +
  • a raster map with maximum water velocity [vmax], + relative flooding intensity [i_vmax], + and the relative time of occurrence [t_vmax],
    • +
  • a raster map with maximum flooding intensity [imax] + and the relative time of occurrence[t_imax],
    • +
  • a raster map with the time of arriving of the Wave-Front + [wavefront]
    • +
+where and the raster maps are coded as "prefix" + "_" + "elapsed seconds": +e.g. mydepth_125.

+ Obviously at least one output map prefix must be specified.
-The unit of measurements of output raster maps are expresssed using the International System (S.I.). -
-
+The unit of measurements of output raster maps are expressed using the +International System (S.I.). - -Options:
-Using a specific flag, the user can obtain another raster map with flow directions that can be visualized using a specific display command -(d.rast.arrow) of the GRASS GIS software.
+

Options

+Using a specific flag, the user can obtain another raster map with flow +directions that can be visualized using a specific display command +(d.rast.arrow) +of the GRASS GIS software.

-Actually two different dam failure type are considered by the command: (i) full breach, (ii) partial breach. +Two different dam failure types are currently considered by the module: +(i) full breach, (ii) partial breach.
- - +Full breach. Partial breach.
-In case of total istantaeous dam break (configuration i), the initial velocity is computed directly applying the SWE at the first time step; -while in case of partial dam breach (configuration ii) the user can choose between don't use any hypothesis, like in the previous configuration, +In case of total instantaneous dam break (configuration i), the initial +velocity is computed directly applying the SWE at the first time step; +while in case of partial dam breach (configuration ii) the user can +choose between don't use any hypothesis, like in the previous configuration, or evaluate the initial velocity using the overflow spillway equation:
- V = 0.4 √(2 g h)
+ V = 0.4 +√ +(2 g h)
where V is the water flow velocity expressed in m/s, g is the gravitational acceleration expressed in m/s2 -and h is the water depth in correspondence of the dam breach expresssed in meters (m).
+and h is the water depth in correspondence of the dam breach +expressed in meters (m).
-Optionally the user may modify the initial timestep used for the numerical solution of the SWE (default value = 0.01 s), nevertheless the timestep [], -and choose a specific failure tipe corresponding to different computational method for the initial velocity estimation. +Optionally the user may modify the initial timestep used for the numerical +solution of the SWE (default value = 0.01 s), nevertheless the +timestep [timestep], +and choose a specific failure type corresponding to different computational +method for the initial velocity estimation.

-
- -
-
-Notes -
-
+

NOTES

(GRASS ANSI C command) -

AUTHORS

-Roberto Marzocchi (e-mail) and Massimiliano Cannata (e-mail). -The GRASS tool was developed by Institute of earth science (IST), -University of applied science of Italian Switzerland (SUPSI), Lugano - Division of geomatics web-page
-Actually the debug is assured by:
- - Gter srl (Genoa, Italy)
- - IST -SUPSI (Lugano, Switzerland)
- - - -The numerical model, originally developed by the National Center for Computational Hydroscience and Engineering of the University of Mississippi, -has been reformulated and modified by the authors introducing important new features to consider the numerical stability and the type of dam failure, -and currently is written in ANSI C programming language within GRASS.

- +

REFERENCES

+[1] Cannata M. & Marzocchi R. (2012). Two-dimensional dam break flooding simulation: a GIS embedded approach. - Natural Hazards 61(3):1143-1159
+[2] Pdf presentation of the work at the +"X Meeting degli Utenti Italiani di GRASS - GFOSS" (It) +web-page
+[3] Pdf presentation of the work at the FOSS4G 2009 (En) - +web-page
+[4] Pdf presentation of the work at the Geoitalia 2011 conference (En) - +document

SEE ALSO

@@ -131,10 +163,35 @@

SEE ALSO

Details of the numerical model are presented in references.
Details of use and developing of are available here.
-

REFERENCES

-[1] Cannata M. & Marzocchi R. (2012). Two-dimensional dam break flooding simulation: a GIS embedded approach. - Natural Hazards 61(3):1143-1159
-[2] Pdf presentation of the work at the "X Meeting degli Utenti Italiani di GRASS - GFOSS" (It) web-page
-[3] Pdf presentation of the work at the FOSS4G 2009 (En) - web-page
-[4] Pdf presentation of the work at the Geoitalia 2011 conference (En)- document
+

AUTHORS

+Roberto Marzocchi (e-mail) +and Massimiliano Cannata (e-mail). +The GRASS tool was developed by Institute of earth science (IST), +University of applied science of Italian Switzerland (SUPSI), Lugano - +Division of geomatics web-page
+Actually the debug is assured by:
+ + + +The numerical model, originally developed by the National Center for +Computational Hydroscience and Engineering of the University of Mississippi, +has been reformulated and modified by the authors introducing important new +features to consider the numerical stability and the type of dam failure, +and currently is written in ANSI C programming language within +GRASS.

+ + \ No newline at end of file From c0e7574c8ed7210d3a820a37d6b5ed10faacbe4e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Edouard=20Choini=C3=A8re?= <27212526+echoix@users.noreply.github.com> Date: Mon, 3 Jan 2022 21:53:51 -0500 Subject: [PATCH 04/17] r.damflood: Translation to English attempt in SWE --- src/raster/r.damflood/SWE.c | 96 +++++++++++++++++++++++++++---------- src/raster/r.damflood/SWE.h | 7 +++ 2 files changed, 79 insertions(+), 24 deletions(-) diff --git a/src/raster/r.damflood/SWE.c b/src/raster/r.damflood/SWE.c index 096b849ce0..cb7eff0cb2 100644 --- a/src/raster/r.damflood/SWE.c +++ b/src/raster/r.damflood/SWE.c @@ -1,3 +1,20 @@ +/**************************************************************************** + * + * MODULE: r.damflood + * AUTHOR: Roberto Marzocchi - roberto.marzocchi[]supsi.ch (2008) + * Massimiliano Cannata - massimiliano.cannata[]supsi.ch (2008) + * PURPOSE: Estimate the area potentially inundated in case of dam breaking + * + * This file handles the Shallow Water Equations + * + * COPYRIGHT: (C) 2008 by Istituto Scienze della Terra-SUPSI + * + * This program is free software under the GNU General Public + * License (>=v2). Read the COPYING file that comes with GRASS + * for details. + * + *****************************************************************************/ + #include #include #include @@ -11,7 +28,7 @@ #include #include -#include "SWE.h" /* specifical dependency to the header file */ +#include "SWE.h" /* Specific dependency to the header file */ @@ -50,7 +67,11 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float /***************************************************/ /* DA METTERE IN UNA ULTERIORE FUNZIONE fall.c */ /* chiamato sia qua che nel main */ - /* controlla Q=0.0 & volume=0.0 */ + /* controlla Q=0.0 & volume=0.0 */ + /* */ + /* TO BE PLACED IN ADDITIONAL FUNCTION IN fall.c */ + /* called both here and in main */ + /* Check (for) Q=0.0 & volume=0.0 */ float Q, vol_res,fall, volume; /***************************************************/ int test; @@ -64,7 +85,7 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float - // DESCRIPTION OF METHOD (italian --> TRASLATE) + // DESCRIPTION OF METHOD (Italian --> TRANSLATE) // primo ciclo: calcolo nuove altezze dell'acqua al tempo t+1: // - a valle della diga applico l'equazione di continuita' delle shallow water // in pratica la nuova altezza e' valutata attraverso un bilancio dei @@ -74,7 +95,18 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float // fisicamente questo porta a una minore realisticita' ma evita le oscillazioni che // sono causa di instabilita' numerica // - nel caso piu' generale si applicano le equazioni a tutto il lago - + // + // (Rough translation): + // First cycle: Calculation of new water heights at time t + 1: + // - Downstream of the dam: Apply continuity equation to shallow water (?) + // In practice, the new height is evaluated through a balance + // of the incoming and outgoing flows in the two main directions + // - Upstream of the dam: + // - In methods 1 and 2: + // - The continuity equation is applied to the volume of the lake + // Physically this leads to a less realistic but avoids + // the oscillations that causes numerical instability. + // - In the more general case the equations are applied to the whole lake for (row = 1; row < nrows-1; row++) { for (col = 1; col < ncols-1; col++) { @@ -133,7 +165,7 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float F = Fdx - Fsx; // dGup =m_v1[row][col] * m_h1[row][col] ;irezione y - // intercella up + // intercell up if (m_v1[row][col]>0 && m_v1[row-1][col]>0) { Gup = m_v1[row][col] * m_h1[row][col]; } else if (m_v1[row][col]<0 && m_v1[row-1][col]<0) { @@ -150,7 +182,7 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float Gup = h_up * v_up; } - // intercella down + // intercell down if (m_v1[row+1][col]>0 && m_v1[row][col]>0) { Gdw = m_v1[row+1][col] * m_h1[row+1][col]; } else if (m_v1[row+1][col]<0 && m_v1[row][col]<0) { @@ -180,11 +212,11 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float } G = Gup - Gdw; - //equazione + //Equation m_h2[row][col] = m_h1[row][col] - timestep / res_ew * F - timestep / res_ns * G; /*if ((row==20||row==21||row==22||row==23)&&(col==18||col==19)){ - printf("EQ. CONTINUITA' --> row:%d, col:%d\n)",row, col); + printf("EQ. CONTINUITY --> row:%d, col:%d\n)",row, col); printf("m_h1[row][col]:%f,m_u1[row][col]:%f,m_v1[row][col]:%f",m_h1[row][col],m_u1[row][col],m_v1[row][col]); printf("m_h1[row][col+1]:%f,m_h1[row][col-1]:%f,m_h1[row+1][col]:%f, m_h1[row-1][col]:%f\n",m_h1[row][col+1],m_h1[row][col-1],m_h1[row+1][col], m_h1[row-1][col]); printf("h_dx:%f, h_sx:%f, h_up%f, h_dw:%f\n",h_dx, h_sx, h_up, h_dw); @@ -207,7 +239,7 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float if (m_h2[row][col]<0){ /*G_warning("At the time %f h is lesser than 0 h(%d,%d)=%f",t, row,col,m_h2[row][col]); printf("row:%d, col:%d, H minore di zero: %.30lf)",row, col, m_h2[row][col]); - printf("DATI:\n"); + printf("DATA:\n"); printf("row:%d,col%d,hmin:%g,h2:%.30lf \n ",row,col,hmin,m_h2[row][col]); printf("m_z[row][col]:%f\n", m_z[row][col]); printf("m_h1[row][col]:%.30lf\n",m_h1[row][col]); @@ -225,12 +257,13 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float m_h2[row][col]=0; } - } // fine continuita' a valle (IF check) + } // fine continuita' a valle (IF check) (end continuity downstream) if (method==1 || method==2){ //******************************************************************* // calcolo portata Q uscente dal lago solo nel caso di Hp stramazzo + // Calculation of flow rate Q coming out of the lake only in the case of Hp weir /* HP: method 1 or 2 */ if (m_DAMBREAK[row][col]>0 ){ if ((m_z[row][col]+m_h1[row][col])>(m_z[row][col+1]+m_h1[row][col+1])){ @@ -256,10 +289,12 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float //***************************************************************************** // abbassamento lago (siccome c'e due volte fare poi una function) + // Lowering of the lake (as there is twice do then a function) //***************************************************************************** if (method==1 || method==2){ /* calcolo l'abbassamento sul lago*/ + /* (Calculation of the lowering of the lake)*/ if (num_cell!=0) { fall = (Q * timestep-vol_res) / (num_cell * res_ew * res_ns); } else { @@ -303,7 +338,7 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float - // DESCRIPTION OF METHOD (italian --> TRASLATE) + // DESCRIPTION OF METHOD (Italian --> TRANSLATE) //********************************************************************************** // terzo ciclo completo sulla matrice: applico le --> // EQUAZIONI DEL MOTO IN DIREZIONE X e Y @@ -312,6 +347,14 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float // NOTA: // u(i,j) e v (i,j) sono le velocita' medie della cella i,j /*******************************************************************/ + //******************************************************************/ + // Third complete cycle over the matrix: Apply --> + // EQUATIONS OF MOTION IN DIRECTIONS X and Y + // and then compute u(t+1) and v(t+1) + // + // NOTE: + // u(i,j) and v(i,j) are the average velocities of cells i,j + /*******************************************************************/ for (row = 1; row < nrows-1; row++) { for (col = 1; col < ncols-1; col++) @@ -320,6 +363,7 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float /**********************************************************************************************************************/ /* EQUAZIONE DEL MOTO IN DIREZIONE X */ + /* (EQUATIONS OF MOTION IN DIRECTION X) */ // right intercell if (m_u1[row][col]>0 && m_u1[row][col+1]>0) { Fdx = m_u1[row][col] * m_u1[row][col] * m_h1[row][col]; @@ -367,7 +411,7 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float F = Fdx - Fsx; //y - // intercella up + // intercell up if (m_v1[row][col]>0 && m_v1[row-1][col]>0) { Gup = m_v1[row][col] * m_u1[row][col] * m_h1[row][col]; } else if (m_v1[row][col]<0 && m_v1[row-1][col]<0) { @@ -385,7 +429,7 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float Gup = h_up * v_up * u_up; } - // intercella down + // intercell down if (m_v1[row+1][col]>0 && m_v1[row][col]>0) { Gdw = m_v1[row+1][col] * m_u1[row+1][col] * m_h1[row+1][col]; } else if (m_v1[row+1][col]<0 && m_v1[row][col]<0) { @@ -417,7 +461,7 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float G = Gup - Gdw; - //courant number --> UPWIND METHOD + //Courant number --> UPWIND METHOD if(m_u1[row][col]>0 && m_u1[row][col+1]>0 && m_u1[row][col-1]>0){ test=1; dZ_dx_down = ( (m_h2[row][col+1] + m_z[row][col+1]) - (m_h2[row][col] + m_z[row][col] )) / res_ew; @@ -485,15 +529,16 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float if (m_DAMBREAK[row][col] > 0){ if ((m_z[row][col]+m_h2[row][col]) > (m_z[row][col+1]+m_h2[row][col+1])) - m_u2[row][col] = velocita_breccia(method,m_h2[row][col]); // velocita' sullo stramazzo + m_u2[row][col] = velocita_breccia(method,m_h2[row][col]); // velocita' sullo stramazzo (velocity on the weir) else if ((m_z[row][col] + m_h2[row][col]) > (m_z[row][col-1] + m_h2[row][col-1])) - m_u2[row][col] = - velocita_breccia(method,m_h2[row][col]); // velocita' sullo stramazzo + m_u2[row][col] = - velocita_breccia(method,m_h2[row][col]); // velocita' sullo stramazzo (velocity on the weir) else m_u2[row][col] = 0.0; }else { m_u2[row][col] = 1.0 / m_h2[row][col] * (m_h1[row][col] * m_u1[row][col] - timestep / res_ew * F - timestep / res_ns * G + timestep * S ); } // no velocita' contro la diga + // (No velocity against the dam) /*if (m_z[row][col+1]> water_elevation && m_u2[row][col]>0) m_u2[row][col]=0.0; if (m_z[row][col-1] > water_elevation && m_u2[row][col]<0) @@ -526,6 +571,7 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float /******************************************************************************************************************************/ /* EQUAZIONE DEL MOTO IN DIREZIONE Y */ + /* (EQUATIONS OF MOTION IN DIRECTION Y) */ // right intercell if (m_u1[row][col]>0 && m_u1[row][col+1]>0) { Fdx = m_u1[row][col] * m_v1[row][col] * m_h1[row][col]; @@ -576,7 +622,7 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float //y - // intercella up + // intercell up if (m_v1[row][col]>0 && m_v1[row-1][col]>0) { Gup = m_v1[row][col] * m_v1[row][col] * m_h1[row][col]; } else if (m_v1[row][col]<0 && m_v1[row-1][col]<0) { @@ -593,7 +639,7 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float Gup = h_up * v_up * v_up; } - // intercella down + // intercell down if (m_v1[row+1][col]>0 && m_v1[row][col]>0) { Gdw = m_v1[row+1][col] * m_v1[row+1][col] * m_h1[row+1][col]; } else if (m_v1[row+1][col]<0 && m_v1[row][col]<0) { @@ -612,7 +658,7 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float if(m_DAMBREAK[row-1][col]>0.0 && ((m_h2[row][col]+m_z[row][col]) < (m_h2[row-1][col]+m_z[row-1][col]))){ - Gup = m_h1[row-1][col]* pow((-velocita_breccia(method,m_h1[row-1][col])),2.0); // -0.4 al quadrato perde il segno meno + Gup = m_h1[row-1][col]* pow((-velocita_breccia(method,m_h1[row-1][col])),2.0); // -0.4 al quadrato perde il segno meno (-0.4 squared loses the minus sign) if(m_h2[row-1][col]==0) Gup=0.0; } @@ -624,7 +670,7 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float G = Gup - Gdw; - //courant number --> UPWIND METHOD + //Courant number --> UPWIND METHOD if (m_v1[row][col]>0 && m_v1[row-1][col]>0 && m_v1[row+1][col]>0){ dZ_dy_down = ((m_h2[row-1][col] + m_z[row-1][col]) - (m_h2[row][col] + m_z[row][col]) ) / res_ns; if (m_h2[row+1][col]==0 && m_z[row+1][col]>(m_h2[row][col] + m_z[row][col])) { @@ -689,9 +735,9 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float if (m_DAMBREAK[row][col] > 0.0 ){ if ((m_z[row][col]+m_h2[row][col]) > (m_z[row-1][col] + m_h2[row-1][col])) - m_v2[row][col] = velocita_breccia(method,m_h2[row][col]); // velocita sullo stramazzo + m_v2[row][col] = velocita_breccia(method,m_h2[row][col]); // velocita sullo stramazzo (velocity on the weir) else if ((m_z[row][col]+m_h2[row][col]) > (m_z[row+1][col] + m_h2[row+1][col])) - m_v2[row][col] = -velocita_breccia(method,m_h2[row][col]); // velocita sullo stramazzo + m_v2[row][col] = -velocita_breccia(method,m_h2[row][col]); // velocita sullo stramazzo (velocity on the weir) else m_v2[row][col] = 0.0; }else{ @@ -699,6 +745,7 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float } // no velocita' contro la diga + // (No velocity against the dam) /*if (m_z[row-1][col] > water_elevation && m_v2[row][col] >0) m_v2[row][col]=0.0; if (m_z[row+1][col] > water_elevation && m_v2[row][col] < 0 ) @@ -726,7 +773,7 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float - //************** stampa ******************************************************** + //************** Prints ******************************************************** //if ((t>6.8 && m_v2[row][col]!=m_v1[row][col]) && (row==87) && (col == 193)) { /*if (fabs(m_v2[row][col])>=1000.0){ G_warning("At the time %f v(%d,%d)=%f", t, row,col,m_v2[row][col]); @@ -736,9 +783,10 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float } else { // tolgo hhmin) + } // ciclo if (h>hmin) (Loop if h>hmin) } } diff --git a/src/raster/r.damflood/SWE.h b/src/raster/r.damflood/SWE.h index c715966808..79e51b4336 100644 --- a/src/raster/r.damflood/SWE.h +++ b/src/raster/r.damflood/SWE.h @@ -6,6 +6,13 @@ originariamente sviluppata per r.damflood (GRASS command) nel caso generico dare una matrice con 2 raster di 0 **m_DAMBREAK & **m_lake e method=3 +returns void +*/ +/* Function to solve Shallow Water Equations +Originally developed for r.damflood (GRASS module) +In the generic case give a matrix with 2 rasters of 0 **m_DAMBREAK & **m_lake +and method=3 + returns void */ From 7809e78ef04f21031432f2373df6dfed274b3c8d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Edouard=20Choini=C3=A8re?= <27212526+echoix@users.noreply.github.com> Date: Sat, 22 Jan 2022 11:55:06 -0500 Subject: [PATCH 05/17] r.damflood: Improve manual and broken link --- src/raster/r.damflood/r.damflood.html | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/raster/r.damflood/r.damflood.html b/src/raster/r.damflood/r.damflood.html index 78c54cc23d..fbad31da7b 100644 --- a/src/raster/r.damflood/r.damflood.html +++ b/src/raster/r.damflood/r.damflood.html @@ -129,8 +129,8 @@

Options

Optionally the user may modify the initial timestep used for the numerical solution of the SWE (default value = 0.01 s), nevertheless the timestep [timestep], -and choose a specific failure type corresponding to different computational -method for the initial velocity estimation. +and choose a specific failure type corresponding to the different computational +methods for the initial velocity estimation.

@@ -156,8 +156,7 @@

REFERENCES

SEE ALSO

r.lake, r.reclass, -d.rast.arrow, -r.inund.fluv. +d.rast.arrow
Details of the numerical model are presented in references.
@@ -193,5 +192,6 @@

AUTHORS

GRASS.

\ No newline at end of file From 283cc2055666ca9c3ba8e5130b1c335dd7362ec2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Edouard=20Choini=C3=A8re?= <27212526+echoix@users.noreply.github.com> Date: Sat, 22 Jan 2022 11:56:00 -0500 Subject: [PATCH 06/17] r.damflood: Translate comments to English --- src/raster/r.damflood/main.c | 119 +++++++++++++++++------------------ 1 file changed, 58 insertions(+), 61 deletions(-) diff --git a/src/raster/r.damflood/main.c b/src/raster/r.damflood/main.c index e94bc64be7..8ccbd5a05d 100644 --- a/src/raster/r.damflood/main.c +++ b/src/raster/r.damflood/main.c @@ -18,7 +18,7 @@ -/* libraries*/ +/* Libraries*/ #include #include #include @@ -32,15 +32,15 @@ #include #include #include -/* function here defined */ -#include "SWE.h" /*function that solves the shallow water equations*/ +/* Function here defined */ +#include "SWE.h" /* Function that solves the shallow water equations*/ //#include -/* simple functions*/ +/* Simple functions*/ #define min(A,B) ((A) < (B) ? (A):(B)) #define max(A,B) ((A) > (B) ? (A):(B)) #define min4(A,B,C,D) min(min(A,B),min(C,D)) @@ -52,8 +52,7 @@ //#define hmin 0.01 -/* lo calcolo dopo in funzione della velocita' massima - (Calculated later as a function of maximum velocity) +/* Calculated later as a function of maximum velocity #define timestep 0.1 */ /* @@ -164,29 +163,29 @@ int main(int argc, char *argv[]){ } startpt; startpt *a_start,*tmp_start;*/ - struct Cell_head cellhd; /* it stores region information and header information of rasters */ + struct Cell_head cellhd; /* Stores region information and header information of rasters */ struct Cell_head window; - //struct History history; /* holds metadata */ + //struct History history; /* Holds metadata */ /* input-output raster file descriptor */ int infd_ELEV,infd_LAKE,infd_DAMBREAK, infd_MANNING, infd_U, infd_V; int outfd_H,outfd_VEL,outfd_VEL_DIR,outfd_HMAX,outfd_T_HMAX,outfd_I_HMAX; int outfd_VMAX,outfd_T_VMAX,outfd_I_VMAX,outfd_DIR_VMAX,outfd_IMAX,outfd_T_IMAX,outfd_WAVEFRONT; - //float g=9.81; define iniziale! + //float g=9.81; Initial definition! - /* mapset name locator */ + /* Mapset name locator */ char *mapset_ELEV,*mapset_LAKE,*mapset_DAMBREAK,*mapset_MANNING,*mapset_U,*mapset_V; - /* buffer for input-output raster */ - FCELL *inrast_ELEV; /* elevation map [m]*/ - DCELL *inrast_LAKE; /* water elevation in the map [m]*/ - FCELL *inrast_DAMBREAK; /* break in the dam*/ - FCELL *inrast_MANNING; /* manning*/ + /* Buffer for input-output raster */ + FCELL *inrast_ELEV; /* Elevation map [m] */ + DCELL *inrast_LAKE; /* Water elevation in the map [m] */ + FCELL *inrast_DAMBREAK; /* Break in the dam */ + FCELL *inrast_MANNING; /* Manning */ DCELL *inrast_U; DCELL *inrast_V; - DCELL *outrast_VEL; /* velocity [m/s] */ + DCELL *outrast_VEL; /* Velocity [m/s] */ DCELL *outrast_VEL_DIR; - DCELL *outrast_H; /* water elevation output [m]*/ + DCELL *outrast_H; /* Water elevation output [m] */ DCELL *outrast_HMAX; DCELL *outrast_I_HMAX; DCELL *outrast_T_HMAX; @@ -197,7 +196,7 @@ int main(int argc, char *argv[]){ DCELL *outrast_IMAX; DCELL *outrast_T_IMAX; DCELL *outrast_WAVEFRONT; - /* cell counters */ + /* Cell counters */ int nrows, ncols; int row, col; int num_cell, num_break; @@ -206,12 +205,12 @@ int main(int argc, char *argv[]){ float Q=0.0, vol_res,fall, volume=0.0; float res_ew ,res_ns; - /* memory matrix */ + /* Memory matrix */ double **m_h1, **m_h2, **m_u1, **m_u2, **m_v1, **m_v2; float **m_z, **m_DAMBREAK, **m_m; float **m_hmax, **m_t_hmax, **m_i_hmax, **m_vmax, **m_t_vmax, **m_i_vmax, **m_dir_vmax, **m_imax, **m_t_imax, **m_wavefront; int ** m_lake; - /* other variables */ + /* Other variables */ double water_elevation, profondita_soglia; double hmin=0.1; int pippo; @@ -231,13 +230,13 @@ int main(int argc, char *argv[]){ int tmp_int; char tmp[15]; /********************/ - // optional instants// + // Optional instants// int ntimes, pp; double *times; double opt_t, time; /********************/ - /* variables to handle user inputs and outputs */ + /* Variables to handle user inputs and outputs */ char *ELEV, *LAKE, *DAMBREAK, *MANNING, *U, *V, *OUT_VEL, *OUT_H, *OUT_HMAX, *OUT_VMAX, *OUT_IMAX, *OUT_WAVEFRONT; @@ -255,7 +254,7 @@ int main(int argc, char *argv[]){ struct { struct Option *met; } opt; - /* initialize GRASS */ + /* Initialize GRASS */ G_gisinit(argv[0]); //pgm_name = argv[0]; @@ -282,10 +281,10 @@ int main(int argc, char *argv[]){ opt.met->answer = "dambreak-without_hypotesis"; //opt.met->guisection = _("Options"); - /* LEGENDA - total_dambreak-without_hypotesis = nessuna ipotesi sulla velocita' iniziale - total_dambreak = Hp altezza critica --> velocita' iniziale (h critica) = 0.93*sqrt(h); - small_dam_breach = Hp stramazzo --> velocita' iniziale = 0.4*sqrt(2*g*h) + /* LEGEND + total_dambreak-without_hypotesis = No hypothesis about initial velocity + total_dambreak = Hp critical height --> initial velocity (critical h) = 0.93*sqrt(h); + small_dam_breach = Hp dam --> initial velocity = 0.4*sqrt(2*g*h) */ input_TIMESTEP = G_define_option(); @@ -427,7 +426,7 @@ int main(int argc, char *argv[]){ /***********************************************************************************************************************/ - /* get entered parameters */ + /* Get entered parameters */ ELEV=input_ELEV->answer; LAKE=input_LAKE->answer; DAMBREAK=input_DAMBREAK->answer; @@ -470,7 +469,7 @@ int main(int argc, char *argv[]){ - /* find maps in mapset */ + /* Find maps in mapset */ mapset_ELEV = (char *) G_find_raster2(ELEV, ""); if (mapset_ELEV == NULL) G_fatal_error (_("cell file [%s] not found"), ELEV); @@ -504,7 +503,7 @@ int main(int argc, char *argv[]){ G_fatal_error (_("cell file [%s] not found"), V); } - /* open input raster files */ + /* Open input raster files */ if ( (infd_ELEV = Rast_open_old (ELEV, mapset_ELEV)) < 0) G_fatal_error (_("Cannot open cell file [%s]"), ELEV); if ( (infd_LAKE = Rast_open_old (LAKE, mapset_LAKE)) < 0) @@ -664,7 +663,7 @@ int main(int argc, char *argv[]){ /* Read every cell in the line buffers */ for (col = 0; col < ncols; col++) { - /* Store values in memory matrix (attenzione valori nulli!) (attention null values!)*/ + /* Store values in memory matrix (attention null values!)*/ m_DAMBREAK[row][col] = ((FCELL *) inrast_DAMBREAK)[col]; m_m[row][col] = ((FCELL *) inrast_MANNING)[col]; m_z[row][col] = ((FCELL *) inrast_ELEV)[col]; @@ -732,7 +731,7 @@ int main(int argc, char *argv[]){ if (method==1 || method==2) { num_cell=0; - /* Search for lakes (cerco il lago) */ + /* Search for lakes */ for (row = 0; row < nrows; row++) { for (col = 0; col < ncols; col++) @@ -751,11 +750,11 @@ int main(int argc, char *argv[]){ { for (col = 0; col < ncols; col++) { - /* Dam break (rottura diga) */ + /* Dam break */ if (m_DAMBREAK[row][col] > 0){ num_break++; G_message("(%d,%d)Cell Dam Breach n° %d",row,col,num_break); - //printf("ho trovato la rottura diga (%d,%d)=\n", row,col); + //printf("I found the dam break (%d,%d)=\n", row,col); //while(!getchar()){ } profondita_soglia = water_elevation-(m_z[row][col] - m_DAMBREAK[row][col]) ; vel_0=velocita_breccia(method,profondita_soglia); @@ -765,15 +764,15 @@ int main(int argc, char *argv[]){ if (vel_0 > vel_max) vel_max=vel_0; } else { - G_fatal_error(_("Didn't find the dambreak. Please select a correct map or adjust the computational region.")); + G_fatal_error(_("Couldn't find the dambreak. Please select a correct map or adjust the computational region.")); } if (m_DAMBREAK[row][col] > 0) { - //cambio il DTM inserendo la rottura della diga + //Change the DTM by inserting the dam break m_z[row][col]=m_z[row][col]-m_DAMBREAK[row][col]; m_lake[row][col]=0; m_h1[row][col]=profondita_soglia; - //printf("la velocità vel_max vale: %f\n",vel_max); + //printf("The velocity vel_max is valid: %f\n",vel_max); //while(!getchar()){ } } } @@ -781,7 +780,6 @@ int main(int argc, char *argv[]){ G_message("The number of lake cell is': %d\n", num_cell); /**************************************/ - /* timestep in funzione di V_0 e res */ /* timestep as a function of V_0 and res */ /**************************************/ //timestep=0.01; @@ -795,7 +793,7 @@ int main(int argc, char *argv[]){ // Uniform drop in of lake (method=1 or method =2) //***************************************************************************** if (method==1 || method==2){ - /* calcolo l'abbassamento sul lago (Calculation of the lowering of the lake) */ + /* Calculation of the lowering of the lake */ if (num_cell!=0) { fall = (volume) / (num_cell * res_ew * res_ns); //printf("volume=%f, fall=%f\n",volume,fall); @@ -807,13 +805,12 @@ int main(int argc, char *argv[]){ for (col = 1; col < ncols-1; col++) { if (m_DAMBREAK[row][col]>0){ - // ragiona se ha senso (I think it makes sense) + // I think it makes sense m_h2[row][col]=m_h1[row][col]-fall; if (m_h2[row][col]<=0) { m_h2[row][col]=0.0; if (m_h1[row][col]>0) { - // questo warning va modificato perchè vale per ogni cella ---> bisogna metterne uno generico che valga quando tutte le celle sono con h=0 - // (This warning must be modified since it is valid for every cell) ---> (You have to put a generic one that is valid when all cells have h=0) + // This warning must be modified since it is valid for every cell ---> You have to put a generic one that is valid when all cells have h=0 num_break--; if (num_break==0){ if (warn1==0){ @@ -832,7 +829,7 @@ int main(int argc, char *argv[]){ num_cell--; } } - }}//end two for cicles + }} //end two for cicles } //end if // There isn't interest to find where is the lake --> everywhere m_lake[row][col]=0 if (method==3){ @@ -840,7 +837,7 @@ int main(int argc, char *argv[]){ { for (col = 0; col < ncols; col++) { - /* Dam break (rottura diga) */ + /* Dam break */ if (m_DAMBREAK[row][col] > 0){ m_z[row][col]=m_z[row][col]-m_DAMBREAK[row][col]; m_DAMBREAK[row][col]=-1.0; @@ -850,7 +847,7 @@ int main(int argc, char *argv[]){ m_lake[row][col]=0; }}}} - G_percent(nrows, nrows, 1); /* finish it */ + G_percent(nrows, nrows, 1); /* Finish it */ G_message("Model running"); @@ -859,7 +856,7 @@ int main(int argc, char *argv[]){ //printf("************************************************\n"); //while(!getchar()){ } //G_percent(t, TSTOP, timestep); - // ciclo sui tempi + // Ciclo sui tempi //G_message("timestep =%f,t=%f",timestep,t); if (t>M*100){ if (M*100!=(m-1)*DELTAT) @@ -875,10 +872,10 @@ int main(int argc, char *argv[]){ shallow_water(m_h1,m_u1,m_v1,m_z,m_DAMBREAK,m_m,m_lake,m_h2,m_u2,m_v2,row,col,nrows,ncols,timestep,res_ew,res_ns,method,num_cell, num_break,t); - //*************************************** overwriting ********************************************* + //*************************************** Overwriting ********************************************* timestep_ct=0; if (t0 || m_v2[nrows-2][col]<0 || m_u1[row][1]<0 || m_u1[row][ncols-2]>0 )) { @@ -889,7 +886,7 @@ int main(int argc, char *argv[]){ } /* velocities at the limit of computational region */ //******************************************************************** - // timestep optimization using the CFL stability condition + // Timestep optimization using the CFL stability condition if (m_h2[row][col]>=hmin ){ timestep_ct_temp = max( (fabs(m_u2[row][col])+sqrt(g*m_h2[row][col]))/res_ew , (fabs(m_v2[row][col])+sqrt(g*m_h2[row][col]))/res_ns ); if(timestep_ct_temp>timestep_ct) { @@ -955,19 +952,19 @@ int main(int argc, char *argv[]){ m_h1[row][col] = m_h2[row][col]; }}} - /*------------------------------ new timestep -------------------------------------*/ + /*------------------------------ New timestep -------------------------------------*/ timestep=0.1/timestep_ct; /*-------------------------------------------------------------------------------------*/ //G_message("timestep =%f,m=%d, t=%f",timestep,m, t); /*----------------------------------------------------------------------------------------------*/ - /* qualcosa non va in questo if */ + /* Something is wrong with this if */ /*----------------------------------------------------------------------------------------------*/ - /* controllo se devo scrivere outputs (Check if we need to write outputs) */ + /* Check if we need to write outputs */ if ((input_DELTAT->answer != NULL)||(parm.opt_t->answer != NULL)) { if ((((m*DELTAT-t) <= timestep && (m*DELTAT) < TSTOP) && (input_DELTAT->answer != NULL)) || ((pp < ntimes && (times[pp]-t) < timestep) && (parm.opt_t->answer != NULL))) { - if ((m*DELTAT-t) <= timestep && m*DELTAT < TSTOP) { /* devo cambiare il nome del raster e aggiungere ogni volta _timestep (We need to change the raster name and add _timestep at each time)*/ + if ((m*DELTAT-t) <= timestep && m*DELTAT < TSTOP) { /* We need to change the raster name and add _timestep at each time */ if (OUT_H) { sprintf(name1,"%s%d",OUT_H,m*DELTAT); } @@ -1044,7 +1041,7 @@ int main(int argc, char *argv[]){ } } /* end_col*/ - /*copia righe (Copy lines)!!! */ + /* Copy lines */ if (OUT_VEL) { Rast_put_d_row(outfd_VEL,outrast_VEL); } @@ -1054,7 +1051,7 @@ int main(int argc, char *argv[]){ if (flag_d->answer) { Rast_put_d_row(outfd_VEL_DIR,outrast_VEL_DIR); } - } /* end row */ + } /* End row */ /* Memory cleanup */ if (OUT_VEL) { G_free(outrast_VEL); @@ -1077,7 +1074,7 @@ int main(int argc, char *argv[]){ } //timestep=0.1/timestep_ct; } - }/* end if write outputs*/ + }/* End if write outputs*/ /* else { //G_message("timestep =%f,t=%f",timestep,t); //pippo=1; @@ -1086,7 +1083,7 @@ int main(int argc, char *argv[]){ } */ //G_message("Function SWE applied - t=%f, TSTOP=%d",t,TSTOP); //G_message("timestep =%f,t=%f",timestep,t); -} /* end time loop*/ +} /* End time loop*/ //******************************************************************* @@ -1184,7 +1181,7 @@ if (OUT_WAVEFRONT){ outrast_WAVEFRONT = Rast_allocate_d_buf(); } -G_percent(nrows, nrows, 1); /* finish it */ +G_percent(nrows, nrows, 1); /* Finish it */ for (row = 0; row < nrows; row++){ G_percent (row, nrows, 2); @@ -1321,7 +1318,7 @@ for (row = 0; row < nrows; row++){ } //G_message("pippo, row=%d e nrows=%d", row, nrows); } // end row -/* chiudi i file (Close files) */ +/* Close files */ if (OUT_H) { @@ -1389,7 +1386,7 @@ if (OUT_WAVEFRONT){ } //************************************************************************ -// da sistemare (TODO/To fix) +// TODO/To fix //************************************************************************ /* Add command line incantation to history file */ //G_short_history(result, "raster", &history); @@ -1397,7 +1394,7 @@ if (OUT_WAVEFRONT){ //G_write_history(result, &history); -/* deallocate memory matrix */ +/* Deallocate memory matrix */ G_free_fmatrix(m_DAMBREAK); G_free_fmatrix(m_m); G_free_fmatrix(m_z); From dfefcbb84dec2cfaef10be504784b174f2b9e678 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Edouard=20Choini=C3=A8re?= <27212526+echoix@users.noreply.github.com> Date: Sat, 22 Jan 2022 12:09:11 -0500 Subject: [PATCH 07/17] r.damflood: Continue translating SWE --- src/raster/r.damflood/SWE.c | 69 ++++++++++--------------------------- src/raster/r.damflood/SWE.h | 7 ---- 2 files changed, 19 insertions(+), 57 deletions(-) diff --git a/src/raster/r.damflood/SWE.c b/src/raster/r.damflood/SWE.c index cb7eff0cb2..414539c394 100644 --- a/src/raster/r.damflood/SWE.c +++ b/src/raster/r.damflood/SWE.c @@ -65,10 +65,6 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float double u_sx, u_dx, v_dx, v_sx, v_up, v_dw, u_up, u_dw; /***************************************************/ - /* DA METTERE IN UNA ULTERIORE FUNZIONE fall.c */ - /* chiamato sia qua che nel main */ - /* controlla Q=0.0 & volume=0.0 */ - /* */ /* TO BE PLACED IN ADDITIONAL FUNCTION IN fall.c */ /* called both here and in main */ /* Check (for) Q=0.0 & volume=0.0 */ @@ -85,18 +81,7 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float - // DESCRIPTION OF METHOD (Italian --> TRANSLATE) - // primo ciclo: calcolo nuove altezze dell'acqua al tempo t+1: - // - a valle della diga applico l'equazione di continuita' delle shallow water - // in pratica la nuova altezza e' valutata attraverso un bilancio dei - // flussi in ingresso e in uscita nelle due direzioni principali - // - a monte delle diga: - // - nel metodo 1 e 2 :l'equazione di continuita' e' applicata al volume del lago - // fisicamente questo porta a una minore realisticita' ma evita le oscillazioni che - // sono causa di instabilita' numerica - // - nel caso piu' generale si applicano le equazioni a tutto il lago - // - // (Rough translation): + // DESCRIPTION OF METHOD // First cycle: Calculation of new water heights at time t + 1: // - Downstream of the dam: Apply continuity equation to shallow water (?) // In practice, the new height is evaluated through a balance @@ -164,7 +149,7 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float } F = Fdx - Fsx; - // dGup =m_v1[row][col] * m_h1[row][col] ;irezione y + // dGup =m_v1[row][col] * m_h1[row][col] ; y direction // intercell up if (m_v1[row][col]>0 && m_v1[row-1][col]>0) { Gup = m_v1[row][col] * m_h1[row][col]; @@ -238,7 +223,7 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float if (m_h2[row][col]<0){ /*G_warning("At the time %f h is lesser than 0 h(%d,%d)=%f",t, row,col,m_h2[row][col]); - printf("row:%d, col:%d, H minore di zero: %.30lf)",row, col, m_h2[row][col]); + printf("row:%d, col:%d, H less than zero: %.30lf)",row, col, m_h2[row][col]); printf("DATA:\n"); printf("row:%d,col%d,hmin:%g,h2:%.30lf \n ",row,col,hmin,m_h2[row][col]); printf("m_z[row][col]:%f\n", m_z[row][col]); @@ -262,7 +247,6 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float if (method==1 || method==2){ //******************************************************************* - // calcolo portata Q uscente dal lago solo nel caso di Hp stramazzo // Calculation of flow rate Q coming out of the lake only in the case of Hp weir /* HP: method 1 or 2 */ if (m_DAMBREAK[row][col]>0 ){ @@ -288,13 +272,11 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float //***************************************************************************** - // abbassamento lago (siccome c'e due volte fare poi una function) // Lowering of the lake (as there is twice do then a function) //***************************************************************************** if (method==1 || method==2){ - /* calcolo l'abbassamento sul lago*/ - /* (Calculation of the lowering of the lake)*/ + /* Calculation of the lowering of the lake*/ if (num_cell!=0) { fall = (Q * timestep-vol_res) / (num_cell * res_ew * res_ns); } else { @@ -338,15 +320,7 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float - // DESCRIPTION OF METHOD (Italian --> TRANSLATE) - //********************************************************************************** - // terzo ciclo completo sulla matrice: applico le --> - // EQUAZIONI DEL MOTO IN DIREZIONE X e Y - // e quindi calcolo u(t+1) e v(t+1) - // - // NOTA: - // u(i,j) e v (i,j) sono le velocita' medie della cella i,j - /*******************************************************************/ + // DESCRIPTION OF METHOD //******************************************************************/ // Third complete cycle over the matrix: Apply --> // EQUATIONS OF MOTION IN DIRECTIONS X and Y @@ -362,8 +336,7 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float if (m_lake[row][col]==0 && m_h2[row][col]>=hmin){ /**********************************************************************************************************************/ - /* EQUAZIONE DEL MOTO IN DIREZIONE X */ - /* (EQUATIONS OF MOTION IN DIRECTION X) */ + /* EQUATIONS OF MOTION IN DIRECTION X */ // right intercell if (m_u1[row][col]>0 && m_u1[row][col+1]>0) { Fdx = m_u1[row][col] * m_u1[row][col] * m_h1[row][col]; @@ -399,7 +372,7 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float } if(m_DAMBREAK[row][col+1]>0 && ((m_h2[row][col]+m_z[row][col]) < (m_h2[row][col+1]+m_z[row][col+1]))){ - Fdx = m_h1[row][col+1]* pow(-velocita_breccia(method,m_h1[row][col+1]),2.0); // -vel al quadrato perde il segno meno + Fdx = m_h1[row][col+1]* pow(-velocita_breccia(method,m_h1[row][col+1]),2.0); // -vel squared looses the negative sign if (m_h2[row][col+1]==0) Fdx=0.0; } @@ -529,16 +502,15 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float if (m_DAMBREAK[row][col] > 0){ if ((m_z[row][col]+m_h2[row][col]) > (m_z[row][col+1]+m_h2[row][col+1])) - m_u2[row][col] = velocita_breccia(method,m_h2[row][col]); // velocita' sullo stramazzo (velocity on the weir) + m_u2[row][col] = velocita_breccia(method,m_h2[row][col]); //velocity on the weir else if ((m_z[row][col] + m_h2[row][col]) > (m_z[row][col-1] + m_h2[row][col-1])) - m_u2[row][col] = - velocita_breccia(method,m_h2[row][col]); // velocita' sullo stramazzo (velocity on the weir) + m_u2[row][col] = - velocita_breccia(method,m_h2[row][col]); //velocity on the weir else m_u2[row][col] = 0.0; }else { m_u2[row][col] = 1.0 / m_h2[row][col] * (m_h1[row][col] * m_u1[row][col] - timestep / res_ew * F - timestep / res_ns * G + timestep * S ); } - // no velocita' contro la diga - // (No velocity against the dam) + // No velocity against the dam /*if (m_z[row][col+1]> water_elevation && m_u2[row][col]>0) m_u2[row][col]=0.0; if (m_z[row][col-1] > water_elevation && m_u2[row][col]<0) @@ -547,7 +519,7 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float if ((timestep/res_ew*(fabs(m_u2[row][col])+sqrt(g*m_h2[row][col])))>1.0){ G_warning("At time %f the Courant-Friedrich-Lewy stability condition isn't respected",t); - /*G_message("velocita' lungo x\n"); + /*G_message("x long velocity \n"); G_message("row:%d, col%d \n",row,col); G_message("dZ_dx_down:%f, dZ_dx_up:%f,cr_up:%f, cr_down:%f\n" , dZ_dx_down,dZ_dx_up, cr_up, cr_down); G_message("Z_piu:%f,Z_meno:%f\n", Z_piu, Z_meno); @@ -570,8 +542,7 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float /******************************************************************************************************************************/ - /* EQUAZIONE DEL MOTO IN DIREZIONE Y */ - /* (EQUATIONS OF MOTION IN DIRECTION Y) */ + /* EQUATIONS OF MOTION IN DIRECTION Y */ // right intercell if (m_u1[row][col]>0 && m_u1[row][col+1]>0) { Fdx = m_u1[row][col] * m_v1[row][col] * m_h1[row][col]; @@ -658,7 +629,7 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float if(m_DAMBREAK[row-1][col]>0.0 && ((m_h2[row][col]+m_z[row][col]) < (m_h2[row-1][col]+m_z[row-1][col]))){ - Gup = m_h1[row-1][col]* pow((-velocita_breccia(method,m_h1[row-1][col])),2.0); // -0.4 al quadrato perde il segno meno (-0.4 squared loses the minus sign) + Gup = m_h1[row-1][col]* pow((-velocita_breccia(method,m_h1[row-1][col])),2.0); // -0.4 squared loses the minus sign if(m_h2[row-1][col]==0) Gup=0.0; } @@ -735,17 +706,16 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float if (m_DAMBREAK[row][col] > 0.0 ){ if ((m_z[row][col]+m_h2[row][col]) > (m_z[row-1][col] + m_h2[row-1][col])) - m_v2[row][col] = velocita_breccia(method,m_h2[row][col]); // velocita sullo stramazzo (velocity on the weir) + m_v2[row][col] = velocita_breccia(method,m_h2[row][col]); // velocity on the weir else if ((m_z[row][col]+m_h2[row][col]) > (m_z[row+1][col] + m_h2[row+1][col])) - m_v2[row][col] = -velocita_breccia(method,m_h2[row][col]); // velocita sullo stramazzo (velocity on the weir) + m_v2[row][col] = -velocita_breccia(method,m_h2[row][col]); // velocity on the weir else m_v2[row][col] = 0.0; }else{ m_v2[row][col] = 1.0 / m_h2[row][col] * (m_h1[row][col] * m_v1[row][col] - timestep / res_ew * F - timestep / res_ns * G + timestep * S); } - // no velocita' contro la diga - // (No velocity against the dam) + // No velocity against the dam /*if (m_z[row-1][col] > water_elevation && m_v2[row][col] >0) m_v2[row][col]=0.0; if (m_z[row+1][col] > water_elevation && m_v2[row][col] < 0 ) @@ -753,7 +723,7 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float if ((timestep/res_ns*(abs(abs(m_v2[row][col])+sqrt(g*m_h2[row][col]))))>1){ G_warning("At time: %f the Courant-Friedrich-Lewy stability condition isn't respected",t); - /*G_message("EQ. MOTO DIR Y' --> row:%d, col:%d\n)",row, col); + /*G_message("EQ. MOTION DIR Y' --> row:%d, col:%d\n)",row, col); G_message("m_h1[row][col]:%f,m_u1[row][col]:%f,m_v1[row][col]:%f",m_h1[row][col],m_u1[row][col],m_v1[row][col]); G_message("m_h1[row][col+1]:%f,m_h1[row][col-1]:%f,m_h1[row+1][col]:%f, m_h1[row-1][col]:%f\n",m_h1[row][col+1],m_h1[row][col-1],m_h1[row+1][col], m_h1[row-1][col]); G_message("h_dx:%f, h_sx:%f, h_up%f, h_dw:%f\n",h_dx, h_sx, h_up, h_dw); @@ -782,11 +752,10 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float } else { - // tolgo hhmin) (Loop if h>hmin) + } // Loop if h>hmin } } diff --git a/src/raster/r.damflood/SWE.h b/src/raster/r.damflood/SWE.h index 79e51b4336..5d78b01570 100644 --- a/src/raster/r.damflood/SWE.h +++ b/src/raster/r.damflood/SWE.h @@ -1,13 +1,6 @@ float velocita_breccia(int i,double h); -/*Funzione per risolvere le shallow water equations -originariamente sviluppata per r.damflood (GRASS command) -nel caso generico dare una matrice con 2 raster di 0 **m_DAMBREAK & **m_lake -e method=3 - -returns void -*/ /* Function to solve Shallow Water Equations Originally developed for r.damflood (GRASS module) In the generic case give a matrix with 2 rasters of 0 **m_DAMBREAK & **m_lake From 3ad438bd8aa98bd1db7fa0871e9d23e1bddaf517 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Edouard=20Choini=C3=A8re?= <27212526+echoix@users.noreply.github.com> Date: Sat, 22 Jan 2022 12:16:41 -0500 Subject: [PATCH 08/17] r.damflood: SWE keep original tab indentation --- src/raster/r.damflood/SWE.c | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/raster/r.damflood/SWE.c b/src/raster/r.damflood/SWE.c index 414539c394..861f5028de 100644 --- a/src/raster/r.damflood/SWE.c +++ b/src/raster/r.damflood/SWE.c @@ -65,8 +65,8 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float double u_sx, u_dx, v_dx, v_sx, v_up, v_dw, u_up, u_dw; /***************************************************/ - /* TO BE PLACED IN ADDITIONAL FUNCTION IN fall.c */ - /* called both here and in main */ + /* TO BE PLACED IN ADDITIONAL FUNCTION IN fall.c */ + /* called both here and in main */ /* Check (for) Q=0.0 & volume=0.0 */ float Q, vol_res,fall, volume; /***************************************************/ @@ -321,14 +321,14 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float // DESCRIPTION OF METHOD - //******************************************************************/ - // Third complete cycle over the matrix: Apply --> - // EQUATIONS OF MOTION IN DIRECTIONS X and Y - // and then compute u(t+1) and v(t+1) - // - // NOTE: - // u(i,j) and v(i,j) are the average velocities of cells i,j - /*******************************************************************/ + //******************************************************************/ + // Third complete cycle over the matrix: Apply --> + // EQUATIONS OF MOTION IN DIRECTIONS X and Y + // and then compute u(t+1) and v(t+1) + // + // NOTE: + // u(i,j) and v(i,j) are the average velocities of cells i,j + /*******************************************************************/ for (row = 1; row < nrows-1; row++) { for (col = 1; col < ncols-1; col++) From 11fea4cd28b7c735400fc7fcec0729d0c4f4ebfd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Edouard=20Choini=C3=A8re?= <27212526+echoix@users.noreply.github.com> Date: Sat, 22 Jan 2022 12:59:43 -0500 Subject: [PATCH 09/17] r.damflood: End r.damflood.html with newline. Co-authored-by: Nicklas Larsson --- src/raster/r.damflood/r.damflood.html | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/raster/r.damflood/r.damflood.html b/src/raster/r.damflood/r.damflood.html index fbad31da7b..9c9cc76b99 100644 --- a/src/raster/r.damflood/r.damflood.html +++ b/src/raster/r.damflood/r.damflood.html @@ -194,4 +194,4 @@

AUTHORS

\ No newline at end of file +--> From f4e304ddb5616c83daa34be545cf93df3828306a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Edouard=20Choini=C3=A8re?= <27212526+echoix@users.noreply.github.com> Date: Wed, 26 Jan 2022 08:51:32 -0500 Subject: [PATCH 10/17] Update src/raster/r.damflood/SWE.c --- src/raster/r.damflood/SWE.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/raster/r.damflood/SWE.c b/src/raster/r.damflood/SWE.c index 861f5028de..44b5a9cc2d 100644 --- a/src/raster/r.damflood/SWE.c +++ b/src/raster/r.damflood/SWE.c @@ -242,7 +242,7 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float m_h2[row][col]=0; } - } // fine continuita' a valle (IF check) (end continuity downstream) + } // end continuity downstream (IF check) if (method==1 || method==2){ From 575f2fc82cf7b85ceaf0493da1e87605a8771fab Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Edouard=20Choini=C3=A8re?= <27212526+echoix@users.noreply.github.com> Date: Sat, 29 Jan 2022 10:35:06 -0500 Subject: [PATCH 11/17] r.damflood: Apply suggestions --- src/raster/r.damflood/SWE.c | 22 +++++++++++----------- src/raster/r.damflood/main.c | 19 +++++++++---------- src/raster/r.damflood/r.damflood.html | 2 +- 3 files changed, 21 insertions(+), 22 deletions(-) diff --git a/src/raster/r.damflood/SWE.c b/src/raster/r.damflood/SWE.c index 44b5a9cc2d..98b0f36360 100644 --- a/src/raster/r.damflood/SWE.c +++ b/src/raster/r.damflood/SWE.c @@ -83,7 +83,7 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float // DESCRIPTION OF METHOD // First cycle: Calculation of new water heights at time t + 1: - // - Downstream of the dam: Apply continuity equation to shallow water (?) + // - Downstream of the dam: Apply continuity equation to shallow water // In practice, the new height is evaluated through a balance // of the incoming and outgoing flows in the two main directions // - Upstream of the dam: @@ -247,8 +247,8 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float if (method==1 || method==2){ //******************************************************************* - // Calculation of flow rate Q coming out of the lake only in the case of Hp weir - /* HP: method 1 or 2 */ + // Calculation of flow rate Q coming out of the lake only in the case of spillway Hp (hypothesis) + /* Hypothesis: method 1 or 2 */ if (m_DAMBREAK[row][col]>0 ){ if ((m_z[row][col]+m_h1[row][col])>(m_z[row][col+1]+m_h1[row][col+1])){ if (t==timestep) @@ -268,15 +268,15 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float } } - }} //end two for cicles (continuity equation) + }} //end two for cycles (continuity equation) //***************************************************************************** - // Lowering of the lake (as there is twice do then a function) + // Lake depth reduction (as there is twice do then a function) //***************************************************************************** if (method==1 || method==2){ - /* Calculation of the lowering of the lake*/ + /* Lake depth reduction calculation*/ if (num_cell!=0) { fall = (Q * timestep-vol_res) / (num_cell * res_ew * res_ns); } else { @@ -502,9 +502,9 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float if (m_DAMBREAK[row][col] > 0){ if ((m_z[row][col]+m_h2[row][col]) > (m_z[row][col+1]+m_h2[row][col+1])) - m_u2[row][col] = velocita_breccia(method,m_h2[row][col]); //velocity on the weir + m_u2[row][col] = velocita_breccia(method,m_h2[row][col]); //velocity on the spillway else if ((m_z[row][col] + m_h2[row][col]) > (m_z[row][col-1] + m_h2[row][col-1])) - m_u2[row][col] = - velocita_breccia(method,m_h2[row][col]); //velocity on the weir + m_u2[row][col] = - velocita_breccia(method,m_h2[row][col]); //velocity on the spillway else m_u2[row][col] = 0.0; }else { @@ -519,7 +519,7 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float if ((timestep/res_ew*(fabs(m_u2[row][col])+sqrt(g*m_h2[row][col])))>1.0){ G_warning("At time %f the Courant-Friedrich-Lewy stability condition isn't respected",t); - /*G_message("x long velocity \n"); + /*G_message("X component of velocity \n"); G_message("row:%d, col%d \n",row,col); G_message("dZ_dx_down:%f, dZ_dx_up:%f,cr_up:%f, cr_down:%f\n" , dZ_dx_down,dZ_dx_up, cr_up, cr_down); G_message("Z_piu:%f,Z_meno:%f\n", Z_piu, Z_meno); @@ -706,9 +706,9 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float if (m_DAMBREAK[row][col] > 0.0 ){ if ((m_z[row][col]+m_h2[row][col]) > (m_z[row-1][col] + m_h2[row-1][col])) - m_v2[row][col] = velocita_breccia(method,m_h2[row][col]); // velocity on the weir + m_v2[row][col] = velocita_breccia(method,m_h2[row][col]); // velocity on the spillway else if ((m_z[row][col]+m_h2[row][col]) > (m_z[row+1][col] + m_h2[row+1][col])) - m_v2[row][col] = -velocita_breccia(method,m_h2[row][col]); // velocity on the weir + m_v2[row][col] = -velocita_breccia(method,m_h2[row][col]); // velocity on the spillway else m_v2[row][col] = 0.0; }else{ diff --git a/src/raster/r.damflood/main.c b/src/raster/r.damflood/main.c index 8ccbd5a05d..4f694cb0d6 100644 --- a/src/raster/r.damflood/main.c +++ b/src/raster/r.damflood/main.c @@ -32,7 +32,7 @@ #include #include #include -/* Function here defined */ +/* Functions defined in here */ #include "SWE.h" /* Function that solves the shallow water equations*/ //#include @@ -171,7 +171,6 @@ int main(int argc, char *argv[]){ int infd_ELEV,infd_LAKE,infd_DAMBREAK, infd_MANNING, infd_U, infd_V; int outfd_H,outfd_VEL,outfd_VEL_DIR,outfd_HMAX,outfd_T_HMAX,outfd_I_HMAX; int outfd_VMAX,outfd_T_VMAX,outfd_I_VMAX,outfd_DIR_VMAX,outfd_IMAX,outfd_T_IMAX,outfd_WAVEFRONT; - //float g=9.81; Initial definition! /* Mapset name locator */ char *mapset_ELEV,*mapset_LAKE,*mapset_DAMBREAK,*mapset_MANNING,*mapset_U,*mapset_V; @@ -283,8 +282,8 @@ int main(int argc, char *argv[]){ /* LEGEND total_dambreak-without_hypotesis = No hypothesis about initial velocity - total_dambreak = Hp critical height --> initial velocity (critical h) = 0.93*sqrt(h); - small_dam_breach = Hp dam --> initial velocity = 0.4*sqrt(2*g*h) + total_dambreak = Hypothesis critical height --> initial velocity (critical h) = 0.93*sqrt(h); + small_dam_breach = Hypothesis spillway --> initial velocity = 0.4*sqrt(2*g*h) */ input_TIMESTEP = G_define_option(); @@ -793,7 +792,7 @@ int main(int argc, char *argv[]){ // Uniform drop in of lake (method=1 or method =2) //***************************************************************************** if (method==1 || method==2){ - /* Calculation of the lowering of the lake */ + /* Lake depth reduction calculation */ if (num_cell!=0) { fall = (volume) / (num_cell * res_ew * res_ns); //printf("volume=%f, fall=%f\n",volume,fall); @@ -829,7 +828,7 @@ int main(int argc, char *argv[]){ num_cell--; } } - }} //end two for cicles + }} //end two for cycles } //end if // There isn't interest to find where is the lake --> everywhere m_lake[row][col]=0 if (method==3){ @@ -856,7 +855,7 @@ int main(int argc, char *argv[]){ //printf("************************************************\n"); //while(!getchar()){ } //G_percent(t, TSTOP, timestep); - // Ciclo sui tempi + // Loop over time //G_message("timestep =%f,t=%f",timestep,t); if (t>M*100){ if (M*100!=(m-1)*DELTAT) @@ -875,14 +874,14 @@ int main(int argc, char *argv[]){ //*************************************** Overwriting ********************************************* timestep_ct=0; if (t0 || m_v2[nrows-2][col]<0 || m_u1[row][1]<0 || m_u1[row][ncols-2]>0 )) { if (reg_lim==0) { G_warning("At the time %.3f the computational region is smaller than inundation",t); reg_lim=1; - } /* Warning message only a time */ + } /* Warning message only a time */ } /* velocities at the limit of computational region */ //******************************************************************** @@ -1386,7 +1385,7 @@ if (OUT_WAVEFRONT){ } //************************************************************************ -// TODO/To fix +// TODO: Fix the following //************************************************************************ /* Add command line incantation to history file */ //G_short_history(result, "raster", &history); diff --git a/src/raster/r.damflood/r.damflood.html b/src/raster/r.damflood/r.damflood.html index 9c9cc76b99..a8fd7b4e7d 100644 --- a/src/raster/r.damflood/r.damflood.html +++ b/src/raster/r.damflood/r.damflood.html @@ -11,7 +11,7 @@

DESCRIPTION

The numerical model solves the conservative form of the shallow water equations (SWE) using a finite volume method (FVM); -the intercell flux is computed by the "upwind method and the water-level +the intercell flux is computed by the upwind method and the water-level gradient is evaluated by weighted average of both upwind and downwind gradient. Additional details of the specific numerical scheme adopted in the model are presented in references [1].
From a5fd00fd73d0c6b9c19d04e17d7d22d0182c388b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Edouard=20Choini=C3=A8re?= <27212526+echoix@users.noreply.github.com> Date: Sat, 30 Dec 2023 16:52:30 +0000 Subject: [PATCH 12/17] Add comments to function parameters in SWE.h --- src/raster/r.damflood/SWE.h | 33 ++++++++++++++++----------------- 1 file changed, 16 insertions(+), 17 deletions(-) diff --git a/src/raster/r.damflood/SWE.h b/src/raster/r.damflood/SWE.h index 9f8fd0b71d..3a25e8bba1 100644 --- a/src/raster/r.damflood/SWE.h +++ b/src/raster/r.damflood/SWE.h @@ -7,21 +7,20 @@ and method=3 returns void */ - void shallow_water( - double **m_h1, double **m_u1, - double **m_v1, /* water depth and velocities of the i step*/ - float **m_z, /* DTM */ - float **m_DAMBREAK, /* DTM changes (e.g. DTM) */ - float **m_m, /* manning coefficient */ - int **m_lake, /* lake filter, default>0, if equal to 0 --> do not apply - swe!!!*/ - double **m_h2, double **m_u2, - double **m_v2, /* water depth and velocities of the i+1 step*/ - int row, int col, int nrows, int ncols, /* matrix size*/ - float timestep, /* timestep (normally optimized with another function) */ - float res_ew, float res_ns, /* grid resolutions*/ - int method, /* default = 3, various hypothesis*/ - int num_cell, int num_break, /* number of cell of lake only in case of - method 1 or 2, elsewhere 0*/ - double t); /* computational instant */ + double **m_h1, /* water depth of the i step */ + double **m_u1, double **m_v1, /* water velocities of the i step */ + float **m_z, /* DTM */ + float **m_DAMBREAK, /* DTM changes (e.g. DTM) */ + float **m_m, /* manning coefficient */ + int **m_lake, // lake filter, default>0, if equal to 0 --> do not apply swe! + double **m_h2, /* water depth of the i+1 step */ + double **m_u2, double **m_v2, /* water velocities of the i+1 step */ + int row, int col, int nrows, int ncols, /* matrix size */ + float timestep, /* timestep (normally optimized with another function) */ + float res_ew, float res_ns, /* grid resolutions */ + int method, /* default = 3, various hypothesis */ + int num_cell, int num_break, /* number of cells of lake only in case of + method 1 or 2, elsewhere 0 */ + double t /* computational instant */ +); From c79c37b004e2273fb7662ce20bdd17b03c3f0b73 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Edouard=20Choini=C3=A8re?= <27212526+echoix@users.noreply.github.com> Date: Sat, 30 Dec 2023 16:53:03 +0000 Subject: [PATCH 13/17] Fix formatting in main.c --- src/raster/r.damflood/main.c | 201 +++++++++++++++++++---------------- 1 file changed, 109 insertions(+), 92 deletions(-) diff --git a/src/raster/r.damflood/main.c b/src/raster/r.damflood/main.c index 191885e5ae..a52cd8812b 100644 --- a/src/raster/r.damflood/main.c +++ b/src/raster/r.damflood/main.c @@ -11,12 +11,12 @@ * COPYRIGHT: (C) 2008 by Istituto Scienze della Terra-SUPSI * * This program is free software under the GNU General Public - * Licence (>=2). Read the file COPYING that cames with GRASS + * Licence (>=2). Read the file COPYING that comes with GRASS * for details. * ***************************************************************************/ -/* Libraries*/ +/* Libraries */ #include #include #include @@ -30,12 +30,13 @@ #include #include #include + /* Functions defined in here */ -#include "SWE.h" /* Function that solves the shallow water equations*/ +#include "SWE.h" /* Function that solves the shallow water equations */ // #include -/* Simple functions*/ +/* Simple functions */ #define min(A, B) ((A) < (B) ? (A) : (B)) #define max(A, B) ((A) > (B) ? (A) : (B)) #define min4(A, B, C, D) min(min(A, B), min(C, D)) @@ -49,8 +50,7 @@ /* Calculated later as a function of maximum velocity #define timestep 0.1 */ -/* - * global function declaration +/* // global function declaration extern CELL f_c(CELL); extern FCELL f_f(FCELL); @@ -69,8 +69,7 @@ FCELL f_calc(FCELL x) DCELL d_calc(DCELL x) { return x; -} - */ +} */ int **G_alloc_imatrix(int rows, int cols) { @@ -143,7 +142,7 @@ void G_free_vector(double *v) return; } -//********************************************************************************************* +//**************************************************************** /* main program */ int main(int argc, char *argv[]) { @@ -190,6 +189,7 @@ int main(int argc, char *argv[]) DCELL *outrast_IMAX; DCELL *outrast_T_IMAX; DCELL *outrast_WAVEFRONT; + /* Cell counters */ int nrows, ncols; int row, col; @@ -205,6 +205,7 @@ int main(int argc, char *argv[]) float **m_hmax, **m_t_hmax, **m_i_hmax, **m_vmax, **m_t_vmax, **m_i_vmax, **m_dir_vmax, **m_imax, **m_t_imax, **m_wavefront; int **m_lake; + /* Other variables */ double water_elevation, profondita_soglia; double hmin = 0.1; @@ -215,29 +216,30 @@ int main(int argc, char *argv[]) int m = 1, M = 1; int i, i_cont; double vel_0 = 0.0, vel_max = 0.0, t; - /**********************************************************************************/ + + /**************************************************************************/ // Parameters for the optimization of timestep using the CFL stability // condition float timestep, velocity; float timestep_ct, timestep_ct_temp; - /**********************************************************************************/ + /**************************************************************************/ + int DELTAT, TSTOP; char name1[20], name2[20], name3[20], name4[20], name5[20], name6[20], name7[20], name8[20], name9[20]; int tmp_int; char tmp[15]; - /********************/ - // Optional instants// + + // Optional instants int ntimes, pp; double *times; double opt_t, time; - /********************/ /* Variables to handle user inputs and outputs */ char *ELEV, *LAKE, *DAMBREAK, *MANNING, *U, *V, *OUT_VEL, *OUT_H, *OUT_HMAX, *OUT_VMAX, *OUT_IMAX, *OUT_WAVEFRONT; - /***********************************************************************************************************************/ + /*************************************************************************/ /* GRASS structure */ struct GModule *module; struct Option *input_ELEV, *input_LAKE, *input_DAMBREAK, *input_MANNING, @@ -251,6 +253,7 @@ int main(int argc, char *argv[]) struct { struct Option *met; } opt; + /* Initialize GRASS */ G_gisinit(argv[0]); @@ -264,7 +267,6 @@ int main(int argc, char *argv[]) _("Estimate the area potentially inundated in case of dam break"); // OPTIONS - flag_d = G_define_flag(); flag_d->key = 'd'; flag_d->description = @@ -284,9 +286,10 @@ int main(int argc, char *argv[]) /* LEGEND total_dambreak-without_hypotesis = No hypothesis about initial velocity - total_dambreak = Hypothesis critical height --> initial velocity (critical h) = - 0.93*sqrt(h); small_dam_breach = Hypothesis spillway --> initial velocity = - 0.4*sqrt(2*g*h) + total_dambreak = Hypothesis critical height + --> initial velocity (critical h) = 0.93*sqrt(h) + small_dam_breach = Hypothesis spillway + --> initial velocity = 0.4*sqrt(2*g*h) */ input_TIMESTEP = G_define_option(); @@ -436,14 +439,14 @@ int main(int argc, char *argv[]) if (G_parser(argc, argv)) exit(EXIT_FAILURE); - /***********************************************************************************************************************/ + /**************************************************************************/ /* Get entered parameters */ ELEV = input_ELEV->answer; LAKE = input_LAKE->answer; DAMBREAK = input_DAMBREAK->answer; MANNING = input_MANNING->answer; sscanf(input_TIMESTEP->answer, "%f", ×tep); - // timestep=input_TIMESTEP->answer; + // timestep = input_TIMESTEP->answer; if (input_U->answer != NULL && input_V->answer != NULL) { U = input_U->answer; V = input_V->answer; @@ -595,11 +598,11 @@ int main(int argc, char *argv[]) inrast_ELEV = Rast_allocate_f_buf(); // inrast_LAKE = Rast_allocate_buf(Rast_map_type(LAKE, mapset_LAKE)); inrast_LAKE = Rast_allocate_d_buf(); - // inrast_DAMBREAK = Rast_allocate_buf(Rast_map_type(DAMBREAK, - // mapset_DAMBREAK)); + // inrast_DAMBREAK = + // Rast_allocate_buf(Rast_map_type(DAMBREAK, mapset_DAMBREAK)); inrast_DAMBREAK = Rast_allocate_f_buf(); - // inrast_MANNING = Rast_allocate_buf(Rast_map_type(MANNING, - // mapset_MANNING)); + // inrast_MANNING = + // Rast_allocate_buf(Rast_map_type(MANNING, mapset_MANNING)); inrast_MANNING = Rast_allocate_f_buf(); if (input_U->answer != NULL && input_V->answer != NULL) { inrast_U = Rast_allocate_d_buf(); @@ -611,6 +614,7 @@ int main(int argc, char *argv[]) else if (input_V->answer != NULL && input_U->answer == NULL) { inrast_V = Rast_allocate_d_buf(); } + /* Get windows rows & cols */ nrows = Rast_window_rows(); ncols = Rast_window_cols(); @@ -668,18 +672,18 @@ int main(int argc, char *argv[]) else if (input_V->answer != NULL && input_U->answer == NULL) { Rast_get_d_row(infd_V, inrast_V, row); } - /*if (G_get_f_raster_row (infd_ELEV, inrast_ELEV, row) < 0) - G_fatal_error (_("Could not read from <%s>"),ELEV); - if (G_get_d_raster_row (infd_LAKE, inrast_LAKE, row) < 0) - G_fatal_error (_("Could not read from <%s>"),LAKE); - if (G_get_f_raster_row (infd_DAMBREAK, inrast_DAMBREAK, row) < 0) - G_fatal_error (_("Could not read from <%s>"),DAMBREAK); - if (G_get_f_raster_row (infd_MANNING, inrast_MANNING, row) < 0) - G_fatal_error (_("Could not read from <%s>"),MANNING);*/ + // if (G_get_f_raster_row(infd_ELEV, inrast_ELEV, row) < 0) + // G_fatal_error(_("Could not read from <%s>"), ELEV); + // if (G_get_d_raster_row(infd_LAKE, inrast_LAKE, row) < 0) + // G_fatal_error(_("Could not read from <%s>"), LAKE); + // if (G_get_f_raster_row(infd_DAMBREAK, inrast_DAMBREAK, row) < 0) + // G_fatal_error(_("Could not read from <%s>"), DAMBREAK); + // if (G_get_f_raster_row(infd_MANNING, inrast_MANNING, row) < 0) + // G_fatal_error(_("Could not read from <%s>"), MANNING); /* Read every cell in the line buffers */ for (col = 0; col < ncols; col++) { - /* Store values in memory matrix (attention null values!)*/ + /* Store values in memory matrix (attention null values!) */ m_DAMBREAK[row][col] = ((FCELL *)inrast_DAMBREAK)[col]; m_m[row][col] = ((FCELL *)inrast_MANNING)[col]; m_z[row][col] = ((FCELL *)inrast_ELEV)[col]; @@ -765,6 +769,7 @@ int main(int argc, char *argv[]) } } num_break = 0; + G_message("Searching dam breach"); for (row = 0; row < nrows; row++) { for (col = 0; col < ncols; col++) { @@ -773,13 +778,14 @@ int main(int argc, char *argv[]) num_break++; G_message("(%d,%d)Cell Dam Breach n° %d", row, col, num_break); - // printf("I found the dam break (%d,%d)=\n", row,col); - // while(!getchar()){ } + // printf("I found the dam break (%d,%d)=\n", row, col); + // while (!getchar()) { + // } profondita_soglia = water_elevation - (m_z[row][col] - m_DAMBREAK[row][col]); vel_0 = velocita_breccia(method, profondita_soglia); volume = volume + profondita_soglia * res_ew * res_ns; - // printf("Q=%.2f\n",Q); + // printf("Q=%.2f\n", Q); // vel_0 = 0.93 * sqrt(profondita_soglia); if (vel_0 > vel_max) vel_max = vel_0; @@ -795,33 +801,35 @@ int main(int argc, char *argv[]) m_z[row][col] = m_z[row][col] - m_DAMBREAK[row][col]; m_lake[row][col] = 0; m_h1[row][col] = profondita_soglia; - // printf("The velocity vel_max is valid: %f\n",vel_max); - // while(!getchar()){ } + // printf("The velocity vel_max is valid: %f\n", vel_max); + // while (!getchar()) { + // } } } } G_message("The number of lake cell is': %d\n", num_cell); - /**************************************/ - /* timestep as a function of V_0 and res */ - /**************************************/ - // timestep=0.01; - // timestep= ((res_ns+res_ew)/2.0)/(vel_max*50.0); + /********************************************* + * timestep as a function of V_0 and res + *********************************************/ + // timestep = 0.01; + // timestep = ((res_ns + res_ew) / 2.0) / (vel_max * 50.0); // DEVELOPEMENT - //***************************************************************************** + //*********************************************************************/ G_message("vel on the dam break cells is %.2f, and timestep is %.2f", vel_max, timestep); } - //***************************************************************************** - // Uniform drop in of lake (method=1 or method =2) - //***************************************************************************** + /*************************************************************************** + * Uniform drop in of lake (method=1 or method=2) + **************************************************************************/ if (method == 1 || method == 2) { /* Lake depth reduction calculation */ if (num_cell != 0) { fall = (volume) / (num_cell * res_ew * res_ns); - // printf("volume=%f, fall=%f\n",volume,fall); - // while(!getchar()){ } + // printf("volume=%f, fall=%f\n", volume, fall); + // while (!getchar()) { + // } } for (row = 1; row < nrows - 1; row++) { @@ -832,8 +840,8 @@ int main(int argc, char *argv[]) if (m_h2[row][col] <= 0) { m_h2[row][col] = 0.0; if (m_h1[row][col] > 0) { - // This warning must be modified since it is valid for every - // cell ---> You have to put a generic one + // This warning must be modified since it is valid + // for every cell ---> You have to put a generic one // that is valid when all cells have h=0 num_break--; if (num_break == 0) { @@ -867,7 +875,7 @@ int main(int argc, char *argv[]) if (m_DAMBREAK[row][col] > 0) { m_z[row][col] = m_z[row][col] - m_DAMBREAK[row][col]; m_DAMBREAK[row][col] = -1.0; - // timestep=0.01; + // timestep = 0.01; m_lake[row][col] = 0; } else { @@ -884,25 +892,25 @@ int main(int argc, char *argv[]) /* Calculate time step loop */ for (t = 0; t <= TSTOP; t += timestep) { // printf("************************************************\n"); - // while(!getchar()){ } + // while (!getchar()) { } // G_percent(t, TSTOP, timestep); - // Loop over time - // G_message("timestep =%f,t=%f",timestep,t); + // Loop over time + // G_message("timestep =%f,t=%f", timestep, t); if (t > M * 100) { if (M * 100 != (m - 1) * DELTAT) G_percent(ceil(t), TSTOP, 2); G_message("t:%d", M * 100); M++; } - // printf("t:%lf\n",t); + // printf("t:%lf\n", t); - // G_message("Function SWE - t=%f, TSTOP=%d",t,TSTOP); + // G_message("Function SWE - t=%f, TSTOP=%d", t, TSTOP); shallow_water(m_h1, m_u1, m_v1, m_z, m_DAMBREAK, m_m, m_lake, m_h2, m_u2, m_v2, row, col, nrows, ncols, timestep, res_ew, res_ns, method, num_cell, num_break, t); - //*************************************** Overwriting ********************************************* + //**************************** Overwriting ****************************/ timestep_ct = 0; if (t < TSTOP) { /* Open new cycle */ @@ -920,7 +928,7 @@ int main(int argc, char *argv[]) } /* Warning message only a time */ } /* velocities at the limit of computational region */ - //******************************************************************** + //*********************************************************/ // Timestep optimization using the CFL stability condition if (m_h2[row][col] >= hmin) { timestep_ct_temp = max( @@ -930,12 +938,15 @@ int main(int argc, char *argv[]) res_ns); if (timestep_ct_temp > timestep_ct) { timestep_ct = timestep_ct_temp; - // G_message("t=%f,row=%d,col=%d,timestep_ct=%f,m_u2=%f - // m_v2=%f - // m_h2=%f",t,row,col,timestep_ct,m_u2[row][col],m_v2[row][col],m_h2[row][col]); + // G_message("t=%f, row=%d, col=%d, " + // "timestep_ct=%f, " + // "m_u2=%f, m_v2=%f, m_h2=%f", + // t, row, col, + // timestep_ct, m_u2[row][col], + // m_v2[row][col], m_h2[row][col]); } } - //******************************************************************** + /**********************************************************/ m_u1[row][col] = m_u2[row][col]; m_v1[row][col] = m_v2[row][col]; @@ -1028,14 +1039,14 @@ int main(int argc, char *argv[]) } } - /*------------------------------ New timestep -------------------------------------*/ + /*--------------------------- New timestep ---------------------------*/ timestep = 0.1 / timestep_ct; - /*-------------------------------------------------------------------------------------*/ - // G_message("timestep =%f,m=%d, t=%f",timestep,m, t); + /*--------------------------------------------------------------------*/ + // G_message("timestep=%f, m=%d, t=%f", timestep, m, t); - /*----------------------------------------------------------------------------------------------*/ - /* Something is wrong with this if */ - /*----------------------------------------------------------------------------------------------*/ + /*-------------------------------------------------------------------- + * Something is wrong with this if + /*--------------------------------------------------------------------*/ /* Check if we need to write outputs */ if ((input_DELTAT->answer != NULL) || (parm.opt_t->answer != NULL)) { @@ -1043,16 +1054,17 @@ int main(int argc, char *argv[]) (input_DELTAT->answer != NULL)) || ((pp < ntimes && (times[pp] - t) < timestep) && (parm.opt_t->answer != NULL))) { - if ((m * DELTAT - t) <= timestep && - m * DELTAT < TSTOP) { /* We need to change the raster name and add _timestep at each time */ + if ((m * DELTAT - t) <= timestep && m * DELTAT < TSTOP) { + /* We need to change the raster name and add _timestep at + * each time */ if (OUT_H) { - sprintf(name1, "%s%d", OUT_H, m*DELTAT); + sprintf(name1, "%s%d", OUT_H, m * DELTAT); } - if (OUT_VEL){ - sprintf(name2,"%s%d",OUT_VEL,m*DELTAT); - sprintf(name3,"%s%s%d","dir_",OUT_VEL,m*DELTAT); + if (OUT_VEL) { + sprintf(name2, "%s%d", OUT_VEL, m * DELTAT); + sprintf(name3, "%s%s%d", "dir_", OUT_VEL, m * DELTAT); } - G_message("Time: %d, writing the output maps",m*DELTAT); + G_message("Time: %d, writing the output maps", m * DELTAT); m++; } else { @@ -1063,6 +1075,7 @@ int main(int argc, char *argv[]) G_message("Time: %lf, writing an optional output maps %d", t, pp); } + /* Controlling if we can write the raster */ if (OUT_H) { if ((outfd_H = Rast_open_new(name1, DCELL_TYPE)) < 0) { @@ -1147,7 +1160,8 @@ int main(int argc, char *argv[]) ((DCELL *)outrast_H)[col] = m_h1[row][col]; } - } /* end_col*/ + } /* end_col */ + /* Copy lines */ if (OUT_VEL) { Rast_put_d_row(outfd_VEL, outrast_VEL); @@ -1159,6 +1173,7 @@ int main(int argc, char *argv[]) Rast_put_d_row(outfd_VEL_DIR, outrast_VEL_DIR); } } /* End row */ + /* Memory cleanup */ if (OUT_VEL) { G_free(outrast_VEL); @@ -1169,6 +1184,7 @@ int main(int argc, char *argv[]) if (flag_d->answer) { G_free(outrast_VEL_DIR); } + /* Close the raster maps */ if (OUT_VEL) { Rast_close(outfd_VEL); @@ -1179,17 +1195,17 @@ int main(int argc, char *argv[]) if (flag_d->answer) { Rast_close(outfd_VEL_DIR); } - // timestep=0.1/timestep_ct; + // timestep = 0.1 / timestep_ct; } - } /* End if write outputs*/ - /* else { - //G_message("timestep =%f,t=%f",timestep,t); - //pippo=1; - G_message("Function SWE - t=%f, TSTOP=%d",t,TSTOP); - //G_percent(t, TSTOP, timestep); - } */ - // G_message("Function SWE applied - t=%f, TSTOP=%d",t,TSTOP); - // G_message("timestep =%f,t=%f",timestep,t); + } /* End if write outputs */ + /* else { + // G_message("timestep=%f, t=%f", timestep, t); + // pippo = 1; + G_message("Function SWE - t=%f, TSTOP=%d", t, TSTOP); + // G_percent(t, TSTOP, timestep); + } */ + // G_message("Function SWE applied - t=%f, TSTOP=%d", t, TSTOP); + // G_message("timestep=%f, t=%f", timestep, t); } /* End time loop*/ //******************************************************************* @@ -1415,12 +1431,13 @@ int main(int argc, char *argv[]) ((DCELL *)outrast_WAVEFRONT)[col] = m_wavefront[row][col]; } } + } /* end_col */ - } /* end_col*/ if (OUT_H) { Rast_put_row(outfd_H, outrast_H, TYPE_DOUBLE); - /*if (G_put_raster_row (outfd_H, outrast_H, TYPE_DOUBLE) < 0) - G_fatal_error (_("Cannot write to <%s>"),name2);*/ + // if (G_put_raster_row(outfd_H, outrast_H, TYPE_DOUBLE) < 0) { + // G_fatal_error(_("Cannot write to <%s>"), name2); + // } } if (OUT_VEL) { Rast_put_row(outfd_VEL, outrast_VEL, TYPE_DOUBLE); @@ -1454,8 +1471,8 @@ int main(int argc, char *argv[]) } // G_message("pippo, row=%d e nrows=%d", row, nrows); } // end row - /* Close files */ + /* Close files */ if (OUT_H) { G_message("Writing the output final map %s", name1); } From d13187d0f36c9d94a2d4f780b1d9ad837dd81b34 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Edouard=20Choini=C3=A8re?= <27212526+echoix@users.noreply.github.com> Date: Sat, 30 Dec 2023 16:54:07 +0000 Subject: [PATCH 14/17] Fix formatting of comments in SWE.c --- src/raster/r.damflood/SWE.c | 358 ++++++++++++++++++++++-------------- 1 file changed, 217 insertions(+), 141 deletions(-) diff --git a/src/raster/r.damflood/SWE.c b/src/raster/r.damflood/SWE.c index 757ab74437..55c6658837 100644 --- a/src/raster/r.damflood/SWE.c +++ b/src/raster/r.damflood/SWE.c @@ -43,7 +43,7 @@ float velocita_breccia(int i, double h) { // double h; // int i; - // float g=9.81; + // float g = 9.81; float v; if (i == 1) { @@ -62,7 +62,7 @@ void shallow_water(double **m_h1, double **m_u1, double **m_v1, float **m_z, int method, int num_cell, int num_break, double t) { - /*FUNCTION VARIABLES*/ + /* FUNCTION VARIABLES */ double h_dx, h_sx, h_up, h_dw, Fup, Fdw, Fdx, Fsx, Gup, Gdw, Gdx, Gsx; double u_sx, u_dx, v_dx, v_sx, v_up, v_dw, u_up, u_dw; @@ -82,23 +82,24 @@ void shallow_water(double **m_h1, double **m_u1, double **m_v1, float **m_z, // DESCRIPTION OF METHOD // First cycle: Calculation of new water heights at time t + 1: - // - Downstream of the dam: Apply continuity equation to shallow water - // In practice, the new height is evaluated through a balance - // of the incoming and outgoing flows in the two main directions - // - Upstream of the dam: - // - In methods 1 and 2: - // - The continuity equation is applied to the volume of the lake - // Physically this leads to a less realistic but avoids - // the oscillations that causes numerical instability. - // - In the more general case the equations are applied to the whole lake + // - Downstream of the dam: Apply continuity equation to shallow water. + // In practice, the new height is evaluated through a balance + // of the incoming and outgoing flows in the two main directions + // - Upstream of the dam: + // - In methods 1 and 2: + // - The continuity equation is applied to the volume of the lake. + // Physically this leads to a less realistic but avoids the + // oscillations that causes numerical instability. + // - In the more general case, the equations are applied to the whole + // lake for (row = 1; row < nrows - 1; row++) { for (col = 1; col < ncols - 1; col++) { if (m_lake[row][col] == 0 && m_DAMBREAK[row][col] <= 0) { - //*******************************************/ - /* CONTINUITY EQUATION --> h(t+1) */ - //*******************************************/ + /******************************************* + * CONTINUITY EQUATION --> h(t+1) + *******************************************/ // x direction // right intercell if (m_u1[row][col] > 0 && m_u1[row][col + 1] > 0) { @@ -168,8 +169,9 @@ void shallow_water(double **m_h1, double **m_u1, double **m_v1, float **m_z, } F = Fdx - Fsx; - // dGup =m_v1[row][col] * m_h1[row][col]; y direction - // intercell up + /* dGup = m_v1[row][col] * m_h1[row][col]; */ + // y direction + // intercell up if (m_v1[row][col] > 0 && m_v1[row - 1][col] > 0) { Gup = m_v1[row][col] * m_h1[row][col]; } @@ -242,67 +244,103 @@ void shallow_water(double **m_h1, double **m_u1, double **m_v1, float **m_z, m_h2[row][col] = m_h1[row][col] - timestep / res_ew * F - timestep / res_ns * G; - /*if ((row==20||row==21||row==22||row==23)&&(col==18||col==19)){ - printf("EQ. CONTINUITY --> row:%d, col:%d\n)",row, - col); - printf("m_h1[row][col]:%f,m_u1[row][col]:%f,m_v1[row][col]:%f",m_h1[row][col],m_u1[row][col],m_v1[row][col]); - printf("m_h1[row][col+1]:%f,m_h1[row][col-1]:%f,m_h1[row+1][col]:%f, - m_h1[row-1][col]:%f\n",m_h1[row][col+1],m_h1[row][col-1],m_h1[row+1][col], - m_h1[row-1][col]); printf("h_dx:%f, h_sx:%f, h_up%f, - h_dw:%f\n",h_dx, h_sx, h_up, h_dw); - printf("m_u1[row][col+1]:%f,m_u1[row][col-1]:%f,m_v1[row+1][col]:%f, - m_v1[row-1][col]:%f\n",m_u1[row][col+1],m_u1[row][col-1],m_v1[row+1][col], - m_v1[row-1][col]); printf("v_up: %f,v_dw:%f,u_dx:%f,u_sx:%f - \n",v_up, v_dw, u_dx, u_sx); printf("Fdx: %f,Fsx: %f, F: %f, - Gup:%f, Gdw:%f, G: %f\n",Fdx, Fsx, F,Gup, Gdw, G); - printf("m_h2(row,col): %f\n \n", m_h2[row][col]); - }*/ + /* if ((row == 20 || row == 21 || row == 22 || row == 23) && + (col == 18 || col == 19)) { + printf("EQ. CONTINUITY --> row:%d, col:%d\n)", row, col); + printf("m_h1[row][col]:%f, " + "m_u1[row][col]:%f, " + "m_v1[row][col]:%f", + m_h1[row][col], m_u1[row][col], m_v1[row][col]); + printf("m_h1[row][col+1]:%f, " + "m_h1[row][col-1]:%f, " + "m_h1[row+1][col]:%f, " + "m_h1[row-1][col]:% f\n", + m_h1[row][col + 1], m_h1[row][col - 1], + m_h1[row + 1][col], m_h1[row - 1][col]); + printf("h_dx:%f, h_sx:%f, h_up%f, h_dw:%f\n", h_dx, h_sx, + h_up, h_dw); + printf("m_u1[row][col+1]:%f, " + "m_u1[row][col-1]:%f, " + "m_v1[row+1][col]:%f, " + "m_v1[row-1][col]:%f\n", + m_u1[row][col + 1], m_u1[row][col - 1], + m_v1[row + 1][col], m_v1[row - 1][col]); + printf("v_up: %f, v_dw:%f, " + "u_dx:%f, u_sx:%f\n", + v_up, v_dw, u_dx, u_sx); + printf("Fdx: %f,Fsx: %f, F: %f, " + "Gup:%f, Gdw:%f, G: %f\n", + Fdx, Fsx, F, Gup, Gdw, G); + printf("m_h2(row,col): %f\n \n", m_h2[row][col]); + } */ - /*if( (row==1 || row==(nrows-2) || col==1 || col==(ncols-2)) && - (m_v2[1][col]>0 || m_v2[nrows-2][col]<0 || m_u1[row][1]<0 || - m_u1[row][ncols-2]>0 )){ if (warn==0){ G_warning("At the time %.3f the - computational region is smaller than inundation",t); warn=1; - G_message("warn=%d",warn); - } - }*/ + /* if (((row == 1 || row == (nrows - 2)) || + (col == 1 || col == (ncols - 2))) && + (m_v2[1][col] > 0 || m_v2[nrows - 2][col] < 0 || + m_u1[row][1] < 0 || m_u1[row][ncols - 2] > 0)) { + if (warn == 0) { + G_warning("At the time %.3f the computational region " + "is smaller than inundation", + t); + warn = 1; + G_message("warn=%d", warn); + } + } */ if (m_h2[row][col] < 0) { - /*G_warning("At the time %f h is lesser than 0 - h(%d,%d)=%f",t, row,col,m_h2[row][col]); printf("row:%d, col:%d, - H less than zero: %.30lf)",row, col, m_h2[row][col]); - printf("DATA:\n"); - printf("row:%d,col%d,hmin:%g,h2:%.30lf \n - ",row,col,hmin,m_h2[row][col]); printf("m_z[row][col]:%f\n", - m_z[row][col]); printf("m_h1[row][col]:%.30lf\n",m_h1[row][col]); - printf("m_u1[row][col]:%.30lf,m_v1[row][col]:%.30lf\n",m_u1[row][col], - m_v1[row][col]); - printf("m_z[row][col+1]:%f,m_z[row][col-1]:%f,m_z[row+1][col]:%f, - m_z[row-1][col]:%f\n",m_z[row][col+1],m_z[row][col-1],m_z[row+1][col], - m_z[row-1][col]); - printf("m_h1[row][col+1]:%.30lf,m_h1[row][col-1]:%.30lf,m_h1[row+1][col]:%.30lf, - m_h1[row-1][col]:%.30lf\n",m_h1[row][col+1],m_h1[row][col-1],m_h1[row+1][col], - m_h1[row-1][col]); printf("h_dx:%f, h_sx:%f, h_up%f, - h_dw:%f\n",h_dx, h_sx, h_up, h_dw); - printf("m_u1[row][col+1]:%.30lf,m_u1[row][col-1]:%.30lf,m_v1[row+1][col]:%.30lf, - m_v1[row-1][col]:%.30lf\n",m_u1[row][col+1],m_u1[row][col-1],m_v1[row+1][col], - m_v1[row-1][col]); printf("timestep: %.30lf, res_ew: %.30lf, - res_ns:%.30lf\n",timestep, res_ew, res_ns); printf("v_up: - %f,v_dw:%f,u_dx:%f,u_sx:%f \n",v_up, v_dw, u_dx, u_sx); - printf("Fdx: %.30lf,Fsx: %.30lf, F: %.30lf, Gup:%.30lf, - Gdw:%.30lf, G: %.30lf\n",Fdx, Fsx, F,Gup, Gdw, G); printf("row: - %d, col %d, m_h1(row,col): %.30lf, m_h2(row,col): %.30lf \n",row, - col,m_h1[row][col], m_h2[row][col]); - printf("m_DAMBREAk(ROW,COL):%f \n",m_DAMBREAK[row][col]); - while(!getchar()){ }*/ + /* G_warning("At the time %f h is lesser than 0 + h(%d,%d)=%f", t, row, col, m_h2[row][col]); printf("row:%d, + col:%d, H less than zero: %.30lf)", row, col, + m_h2[row][col]); printf("DATA:\n"); + printf("row:%d,col%d,hmin:%g,h2:%.30lf\n", row, col, hmin, + m_h2[row][col]); + printf("m_z[row][col]:%f\n", m_z[row][col]); + printf("m_h1[row][col]:%.30lf\n", m_h1[row][col]); + printf("m_u1[row][col]:%.30lf, m_v1[row][col]:%.30lf\n", + m_u1[row][col], m_v1[row][col]); + printf("m_z[row][col+1]:%f, " + "m_z[row][col-1]:%f, " + "m_z[row+1][col]:%f, " + "m_z[row-1][col]:%f\n", + m_z[row][col + 1], m_z[row][col - 1], + m_z[row + 1][col], m_z[row - 1][col]); + printf("m_h1[row][col+1]:%.30lf, " + "m_h1[row][col-1]:%.30lf, " + "m_h1[row+1][col]:%.30lf, " + "m_h1[row-1][col]:%.30lf\n", + m_h1[row][col + 1], m_h1[row][col - 1], + m_h1[row + 1][col], m_h1[row - 1][col]); + printf("h_dx:%f, h_sx:%f, " + "h_up%f, h_dw:%f\n", + h_dx, h_sx, h_up, h_dw); + printf("m_u1[row][col+1]:%.30lf, " + "m_u1[row][col-1]:%.30lf, " + "m_v1[row+1][col]:%.30lf, " + "m_v1[row-1][col]:%.30lf\n", + m_u1[row][col + 1], m_u1[row][col - 1], + m_v1[row + 1][col], m_v1[row - 1][col]); + printf("timestep: %.30lf, res_ew: %.30lf, res_ns:%.30lf\n", + timestep, res_ew, res_ns); + printf("v_up: %f, v_dw:%f, u_dx:%f, u_sx:%f \n", v_up, v_dw, + u_dx, u_sx); + printf("Fdx: %.30lf, Fsx: %.30lf, F: %.30lf, " + "Gup:%.30lf, Gdw:%.30lf, G: %.30lf\n", + Fdx, Fsx, F, Gup, Gdw, G); + printf("row: %d, col %d, " + "m_h1(row,col): %.30lf, m_h2(row,col): %.30lf \n", + row, col, m_h1[row][col], m_h2[row][col]); + printf("m_DAMBREAK(ROW,COL):%f \n", m_DAMBREAK[row][col]); + while (!getchar()) { + } */ m_h2[row][col] = 0; } } // end continuity downstream (IF check) if (method == 1 || method == 2) { - //******************************************************************* - // Calculation of flow rate Q coming out of the lake only in the case of spillway Hp (hypothesis) - /* Hypothesis: method 1 or 2 */ + // Calculation of flow rate Q coming out of the lake only in the + // case of spillway Hp (hypothesis) + /* Hypothesis: method 1 or 2 */ if (m_DAMBREAK[row][col] > 0) { if ((m_z[row][col] + m_h1[row][col]) > (m_z[row][col + 1] + m_h1[row][col + 1])) { @@ -353,10 +391,10 @@ void shallow_water(double **m_h1, double **m_u1, double **m_v1, float **m_z, fall = (Q * timestep - vol_res) / (num_cell * res_ew * res_ns); } else { - // if (warn1==0){ + // if (warn1 == 0) { G_warning("At the time %.0fs no water go out from lake", t); - // warn1=1; - //} + // warn1 = 1; + // } } vol_res = 0.0; Q = 0.0; @@ -369,10 +407,10 @@ void shallow_water(double **m_h1, double **m_u1, double **m_v1, float **m_z, m_h2[row][col] = 0.0; if (m_h1[row][col] > 0) { num_break--; - /*if (num_break==0){ - G_warning("At the time %.0fs no water go out - from lake",t); - }*/ + /* if (num_break == 0) { + G_warning("At the time %.0fs no water go out " + "from lake", t); + } */ } } } @@ -403,7 +441,7 @@ void shallow_water(double **m_h1, double **m_u1, double **m_v1, float **m_z, for (col = 1; col < ncols - 1; col++) { if (m_lake[row][col] == 0 && m_h2[row][col] >= hmin) { - /**********************************************************************************************************************/ + /*******************************************************/ /* EQUATIONS OF MOTION IN DIRECTION X */ // right intercell if (m_u1[row][col] > 0 && m_u1[row][col + 1] > 0) { @@ -552,7 +590,7 @@ void shallow_water(double **m_h1, double **m_u1, double **m_v1, float **m_z, } G = Gup - Gdw; - // Courant number --> UPWIND METHOD + // Courant number --> UPWIND METHOD if (m_u1[row][col] > 0 && m_u1[row][col + 1] > 0 && m_u1[row][col - 1] > 0) { test = 1; @@ -670,12 +708,12 @@ void shallow_water(double **m_h1, double **m_u1, double **m_v1, float **m_z, (m_z[row][col + 1] + m_h2[row][col + 1])) m_u2[row][col] = velocita_breccia( method, - m_h2[row][col]); //velocity on the spillway + m_h2[row][col]); // velocity on the spillway else if ((m_z[row][col] + m_h2[row][col]) > (m_z[row][col - 1] + m_h2[row][col - 1])) m_u2[row][col] = -velocita_breccia( method, - m_h2[row][col]); //velocity on the spillway + m_h2[row][col]); // velocity on the spillway else m_u2[row][col] = 0.0; } @@ -686,48 +724,62 @@ void shallow_water(double **m_h1, double **m_u1, double **m_v1, float **m_z, timestep / res_ns * G + timestep * S); } // No velocity against the dam - /*if (m_z[row][col+1]> water_elevation && m_u2[row][col]>0) - m_u2[row][col]=0.0; - if (m_z[row][col-1] > water_elevation && m_u2[row][col]<0) - m_u2[row][col]=0.0;*/ + // if (m_z[row][col + 1] > water_elevation + // && m_u2[row][col] > 0) { + // m_u2[row][col] = 0.0; + // } + // if (m_z[row][col - 1] > water_elevation + // && m_u2[row][col] < 0) { + // m_u2[row][col] = 0.0; + // } if ((timestep / res_ew * (fabs(m_u2[row][col]) + sqrt(g * m_h2[row][col]))) > 1.0) { G_warning("At time %f the Courant-Friedrich-Lewy stability " "condition isn't respected", t); - /*G_message("X component of velocity \n"); - G_message("row:%d, col%d \n",row,col); - G_message("dZ_dx_down:%f, dZ_dx_up:%f,cr_up:%f, - cr_down:%f\n" , dZ_dx_down,dZ_dx_up, cr_up, cr_down); - G_message("Z_piu:%f,Z_meno:%f\n", Z_piu, Z_meno); - G_message("dZ_dx:%f\n",dZ_dx); - G_message("m_h1[row][col]:%f, m_h2[row][col]:%f, - m_z[row][col]:%f\n",m_h1[row][col], m_h2[row][col], - m_z[row][col]); G_message("m_h2[row][col+1]:%f, - m_z[row][col+1]:%f,m_h2[row][col-1]:%f, m_z[row][col-1]:%f - \n",m_h2[row][col+1], - m_z[row][col+1],m_h2[row][col-1],m_z[row][col-1]); - G_message("m_h2[row][col+1]:%f,m_h2[row][col-1]:%f,\n",m_h2[row][col+1],m_h2[row][col-1]); - G_message("h_up: %f,h_dw:%f,h_dx:%f,h_sx:%f \n",h_up, h_dw, - h_dx, h_sx); G_message("Fdx: %f,Fsx: %f, F: %f, Gup:%f, - Gdw:%f, G: %.60lf, S: %.60lf \n",Fdx, Fsx, F,Gup, Gdw, G, - S); G_message("m_u1[row][col-1]:%f, m_h1[row][col-1]:%f, - m_u1[row][col+1]:%f, - m_h1[row][col+1]:%f\n",m_u1[row][col-1], m_h1[row][col-1], - m_u1[row][col+1], m_h1[row][col+1]); G_message("timestep:%f, - res_ew:%f\n",timestep,res_ew); - G_message("m_u2[row][col]:%f,m_u1[row][col]:%f\n\n", - m_u2[row][col],m_u1[row][col]); G_warning(" ");*/ + /* G_message("X component of velocity \n"); + G_message("row:%d, col%d \n", row, col); + G_message( + "dZ_dx_down:%f, dZ_dx_up:%f, cr_up:%f, cr_down:%f\n", + dZ_dx_down, dZ_dx_up, cr_up, cr_down); + G_message("Z_piu:%f, Z_meno:%f\n", Z_piu, Z_meno); + G_message("dZ_dx:%f\n", dZ_dx); + G_message("m_h1[row][col]:%f, " + "m_h2[row][col]:%f, " + "m_z[row][col]:%f\n", + m_h1[row][col], m_h2[row][col], m_z[row][col]); + G_message("m_h2[row][col+1]:%f, " + "m_z[row][col+1]:%f, " + "m_h2[row][col-1]:%f, " + "m_z[row][col-1]:%f\n", + m_h2[row][col + 1], m_z[row][col + 1], + m_h2[row][col - 1], m_z[row][col - 1]); + G_message("m_h2[row][col+1]:%f, " + "m_h2[row][col-1]:%f,\n", + m_h2[row][col + 1], m_h2[row][col - 1]); + G_message("h_up: %f, h_dw:%f, h_dx:%f, h_sx:%f\n", h_up, + h_dw, h_dx, h_sx); + G_message("Fdx: %f,Fsx: %f, F: %f, " + "Gup:%f, Gdw:%f, G: %.60lf, S: %.60lf \n", + Fdx, Fsx, F, Gup, Gdw, G, S); + G_message("m_u1[row][col-1]:%f, m_h1[row][col-1]:%f, " + "m_u1[row][col+1]:%f, m_h1[row][col+1]:%f\n", + m_u1[row][col - 1], m_h1[row][col - 1], + m_u1[row][col + 1], m_h1[row][col + 1]); + G_message("timestep:%f, res_ew:%f\n", timestep, res_ew); + G_message("m_u2[row][col]:%f, m_u1[row][col]:%f\n\n", + m_u2[row][col], m_u1[row][col]); + G_warning(" "); */ } if (fabs(m_u2[row][col] >= 1000)) { G_warning("At the time %f u(%d,%d)=%f", t, row, col, m_u2[row][col]); } - /******************************************************************************************************************************/ + /*************************************************************/ - /******************************************************************************************************************************/ + /*************************************************************/ /* EQUATIONS OF MOTION IN DIRECTION Y */ // right intercell if (m_u1[row][col] > 0 && m_u1[row][col + 1] > 0) { @@ -1002,10 +1054,14 @@ void shallow_water(double **m_h1, double **m_u1, double **m_v1, float **m_z, } // No velocity against the dam - /*if (m_z[row-1][col] > water_elevation && m_v2[row][col] >0) - m_v2[row][col]=0.0; - if (m_z[row+1][col] > water_elevation && m_v2[row][col] < 0 ) - m_v2[row][col]=0.0;*/ + // if (m_z[row - 1][col] > water_elevation + // && m_v2[row][col] > 0) { + // m_v2[row][col] = 0.0; + // } + // if (m_z[row + 1][col] > water_elevation + // && m_v2[row][col] < 0) { + // m_v2[row][col] = 0.0; + // } if ((timestep / res_ns * (abs(abs(m_v2[row][col]) + sqrt(g * m_h2[row][col])))) > @@ -1013,38 +1069,58 @@ void shallow_water(double **m_h1, double **m_u1, double **m_v1, float **m_z, G_warning("At time: %f the Courant-Friedrich-Lewy " "stability condition isn't respected", t); - /*G_message("EQ. MOTION DIR Y' --> row:%d, col:%d\n)",row, - col); - G_message("m_h1[row][col]:%f,m_u1[row][col]:%f,m_v1[row][col]:%f",m_h1[row][col],m_u1[row][col],m_v1[row][col]); - G_message("m_h1[row][col+1]:%f,m_h1[row][col-1]:%f,m_h1[row+1][col]:%f, - m_h1[row-1][col]:%f\n",m_h1[row][col+1],m_h1[row][col-1],m_h1[row+1][col], - m_h1[row-1][col]); G_message("h_dx:%f, h_sx:%f, h_up%f, - h_dw:%f\n",h_dx, h_sx, h_up, h_dw); - G_message("m_u1[row][col+1]:%f,m_u1[row][col-1]:%f,m_v1[row+1][col]:%f, - m_v1[row-1][col]:%f\n",m_u1[row][col+1],m_u1[row][col-1],m_v1[row+1][col], - m_v1[row-1][col]); G_message("v_up: - %f,v_dw:%f,u_dx:%f,u_sx:%f \n",v_up, v_dw, u_dx, u_sx); - G_message("Fdx: %f,Fsx: %f, F: %f, Gup:%f, Gdw:%f, G: - %f\n",Fdx, Fsx, F,Gup, Gdw, G); - G_message("m_h2[row][col+1]:%f,m_h2[row][col-1]:%f,m_h2[row+1][col]:%f, - m_h2[row-1][col]:%f\n",m_h1[row][col+1],m_h1[row][col-1],m_h1[row+1][col], - m_h1[row-1][col]); G_message("dZ_dy_down:%f, - dZ_dy_up:%f,cr_up:%f, cr_down:%f\n" , dZ_dy_down,dZ_dy_up, - cr_up, cr_down); G_message("Z_piu:%f,Z_meno:%f\n", Z_piu, - Z_meno); G_message("dZ_dy:%f,\n",dZ_dy); - G_message("u:%f,v:%f,V:%f\n", u,v,V); - G_message("R_i:%f,manning[row][col]:%f\n", R_i, - m_m[row][col]); G_message("S=%f\n",S); + /* G_message("EQ. MOTION DIR Y' --> row:%d, col:%d\n)", row, + col); + G_message("m_h1[row][col]:%f, " + "m_u1[row][col]:%f, " + "m_v1[row][col]:%f", + m_h1[row][col], m_u1[row][col], m_v1[row][col]); + G_message("m_h1[row][col+1]:%f, " + "m_h1[row][col-1]:%f, " + "m_h1[row+1][col]:%f, " + "m_h1[row-1][col]:%f\n", + m_h1[row][col + 1], m_h1[row][col - 1], + m_h1[row + 1][col], m_h1[row - 1][col]); + G_message("h_dx:%f, h_sx:%f, h_up%f, h_dw:%f\n", h_dx, h_sx, + h_up, h_dw); + G_message("m_u1[row][col+1]:%f, " + "m_u1[row][col-1]:%f, " + "m_v1[row+1][col]:%f, " + "m_v1[row-1][col]:%f\n", + m_u1[row][col + 1], m_u1[row][col - 1], + m_v1[row + 1][col], m_v1[row - 1][col]); + G_message("v_up: %f,v_dw:%f,u_dx:%f,u_sx:%f \n", v_up, v_dw, + u_dx, u_sx); + G_message("Fdx: %f,Fsx: %f, F: %f, " + "Gup:%f, Gdw:%f, G: %f\n", + Fdx, Fsx, F, Gup, Gdw, G); + G_message("m_h2[row][col+1]:%f, " + "m_h2[row][col-1]:%f, " + "m_h2[row+1][col]:%f, " + "m_h2[row-1][col]:%f\n", + m_h1[row][col + 1], m_h1[row][col - 1], + m_h1[row + 1][col], m_h1[row - 1][col]); + G_message("dZ_dy_down:%f, dZ_dy_up:%f, " + "cr_up:%f, cr_down:%f\n", + dZ_dy_down, dZ_dy_up, cr_up, cr_down); + G_message("Z_piu:%f, Z_meno:%f\n", Z_piu, Z_meno); + G_message("dZ_dy:%f,\n", dZ_dy); + G_message("u:%f, v:%f, V:%f\n", u, v, V); + G_message("R_i:%f, manning[row][col]:%f\n", R_i, + m_m[row][col]); + G_message("S=%f\n", S); G_message("m_v2(row,col): %f\n \n", m_v2[row][col]); - G_warning(" ");*/ + G_warning(" "); */ } - //************** Prints ******************************************************** - //if ((t>6.8 && m_v2[row][col]!=m_v1[row][col]) && (row==87) && (col == 193)) { - /*if (fabs(m_v2[row][col])>=1000.0){ - G_warning("At the time %f v(%d,%d)=%f", t, - row,col,m_v2[row][col]); - }*/ + /*************************** Prints ***************************/ + // if ((t > 6.8 && m_v2[row][col] != m_v1[row][col]) && + // (row == 87) && (col == 193)) { + /* if (fabs(m_v2[row][col]) >= 1000.0) { + G_warning("At the time %f v(%d,%d)=%f", + t, row, col, m_v2[row][col]); + } */ + // } } else { // Remove h Date: Sat, 30 Dec 2023 18:59:08 +0000 Subject: [PATCH 15/17] r.damflood: Fix typos and formatting of manual --- src/raster/r.damflood/r.damflood.html | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/src/raster/r.damflood/r.damflood.html b/src/raster/r.damflood/r.damflood.html index d5969a7cde..99434810b6 100644 --- a/src/raster/r.damflood/r.damflood.html +++ b/src/raster/r.damflood/r.damflood.html @@ -1,6 +1,5 @@

DESCRIPTION

- r.damflood - The definition of flooding areas is of considerable importance for both the risk analysis and the emergency management. This module, in particular, is an embedded GRASS GIS hydrodynamic 2D model @@ -159,14 +158,14 @@

SEE ALSO


Details of the numerical model are presented in references.
-Details of use and developing of are available here.
+Usage and developping details of are available here.

AUTHORS

Roberto Marzocchi (e-mail) and Massimiliano Cannata (e-mail). -The GRASS tool was developed by Institute of earth science (IST), +The GRASS tool was developed by the Institute of earth science (IST), University of applied science of Italian Switzerland (SUPSI), Lugano - Division of geomatics web-page
@@ -189,8 +188,3 @@

AUTHORS

features to consider the numerical stability and the type of dam failure, and currently is written in ANSI C programming language within GRASS.

- - From db2c0d8f1c5546d1462d31cadbaf499ea975939f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Edouard=20Choini=C3=A8re?= <27212526+echoix@users.noreply.github.com> Date: Sat, 30 Dec 2023 19:49:47 -0500 Subject: [PATCH 16/17] Apply suggestions from code review Co-authored-by: Veronica Andreo --- src/raster/r.damflood/SWE.c | 4 ++-- src/raster/r.damflood/main.c | 10 +++++----- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/raster/r.damflood/SWE.c b/src/raster/r.damflood/SWE.c index 55c6658837..faf1a6b5dd 100644 --- a/src/raster/r.damflood/SWE.c +++ b/src/raster/r.damflood/SWE.c @@ -88,7 +88,7 @@ void shallow_water(double **m_h1, double **m_u1, double **m_v1, float **m_z, // - Upstream of the dam: // - In methods 1 and 2: // - The continuity equation is applied to the volume of the lake. - // Physically this leads to a less realistic but avoids the + // Physically this leads to less realism but avoids the // oscillations that causes numerical instability. // - In the more general case, the equations are applied to the whole // lake @@ -435,7 +435,7 @@ void shallow_water(double **m_h1, double **m_u1, double **m_v1, float **m_z, // and then compute u(t+1) and v(t+1) // // NOTE: - // u(i,j) and v(i,j) are the average velocities of cells i,j + // u(i,j) and v(i,j) are the average velocities of cell i,j /*******************************************************************/ for (row = 1; row < nrows - 1; row++) { for (col = 1; col < ncols - 1; col++) { diff --git a/src/raster/r.damflood/main.c b/src/raster/r.damflood/main.c index a52cd8812b..f3c49e3987 100644 --- a/src/raster/r.damflood/main.c +++ b/src/raster/r.damflood/main.c @@ -792,7 +792,7 @@ int main(int argc, char *argv[]) } else { G_fatal_error( - _("Couldn't find the dambreak. Please select a correct " + _("Couldn't find the dambreak. Please select the correct " "map or adjust the computational region.")); } @@ -835,14 +835,14 @@ int main(int argc, char *argv[]) for (row = 1; row < nrows - 1; row++) { for (col = 1; col < ncols - 1; col++) { if (m_DAMBREAK[row][col] > 0) { - // I think it makes sense + // Checks if it makes sense m_h2[row][col] = m_h1[row][col] - fall; if (m_h2[row][col] <= 0) { m_h2[row][col] = 0.0; if (m_h1[row][col] > 0) { - // This warning must be modified since it is valid - // for every cell ---> You have to put a generic one - // that is valid when all cells have h=0 + // This warning is modified since it is valid + // for every cell ---> A generic one is needed + // for the case when all cells have h=0 num_break--; if (num_break == 0) { if (warn1 == 0) { From 24c42d6be0f9de4b184ef88323bcdd624365a11b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Edouard=20Choini=C3=A8re?= <27212526+echoix@users.noreply.github.com> Date: Sat, 30 Dec 2023 20:07:16 -0500 Subject: [PATCH 17/17] Update main.c --- src/raster/r.damflood/main.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/raster/r.damflood/main.c b/src/raster/r.damflood/main.c index f3c49e3987..993037bc12 100644 --- a/src/raster/r.damflood/main.c +++ b/src/raster/r.damflood/main.c @@ -792,8 +792,8 @@ int main(int argc, char *argv[]) } else { G_fatal_error( - _("Couldn't find the dambreak. Please select the correct " - "map or adjust the computational region.")); + _("Couldn't find the dambreak. Please select the " + "correct map or adjust the computational region.")); } if (m_DAMBREAK[row][col] > 0) {