Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

r.damflood: Update manual and migrate Italian to English #683

Open
wants to merge 26 commits into
base: grass8
Choose a base branch
from
Open

r.damflood: Update manual and migrate Italian to English #683

Show file tree
Hide file tree
Changes from 12 commits
Commits
Show all changes
26 commits
Select commit Hold shift + click to select a range
9ac2a06
Add english comments
echoix Jan 3, 2022
7eade51
Spelling fixes in r.damflood
echoix Jan 3, 2022
d80bfd3
r.damflood manual page revision
echoix Jan 3, 2022
c0e7574
r.damflood: Translation to English attempt in SWE
echoix Jan 4, 2022
b74a1b9
Merge branch 'OSGeo:grass8' into r.damflood-patch-1
echoix Jan 16, 2022
bde6fbf
Merge branch 'OSGeo:grass8' into r.damflood-patch-1
echoix Jan 16, 2022
7809e78
r.damflood: Improve manual and broken link
echoix Jan 22, 2022
283cc20
r.damflood: Translate comments to English
echoix Jan 22, 2022
af848a3
Merge branch 'r.damflood-patch-1' of https://github.com/echoix/grass-…
echoix Jan 22, 2022
dfefcbb
r.damflood: Continue translating SWE
echoix Jan 22, 2022
3ad438b
r.damflood: SWE keep original tab indentation
echoix Jan 22, 2022
11fea4c
r.damflood: End r.damflood.html with newline.
echoix Jan 22, 2022
f4e304d
Update src/raster/r.damflood/SWE.c
echoix Jan 26, 2022
575f2fc
r.damflood: Apply suggestions
echoix Jan 29, 2022
796b43c
Merge branch 'grass8' into r.damflood-patch-1
echoix Dec 30, 2023
e3ca1ef
Merge branch 'grass8' into r.damflood-patch-1
echoix Dec 30, 2023
60d3351
Merge branch 'grass8' into r.damflood-patch-1
echoix Dec 30, 2023
5bd8815
Merge branch 'grass8' into r.damflood-patch-1
echoix Dec 30, 2023
f0ebe9a
Merge branch 'grass8' into r.damflood-patch-1
echoix Dec 30, 2023
b596a54
Merge branch 'grass8' into r.damflood-patch-1
echoix Dec 30, 2023
a5fd00f
Add comments to function parameters in SWE.h
echoix Dec 30, 2023
c79c37b
Fix formatting in main.c
echoix Dec 30, 2023
d13187d
Fix formatting of comments in SWE.c
echoix Dec 30, 2023
76149f0
r.damflood: Fix typos and formatting of manual
echoix Dec 30, 2023
db2c0d8
Apply suggestions from code review
echoix Dec 31, 2023
24c42d6
Update main.c
echoix Dec 31, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
125 changes: 71 additions & 54 deletions src/raster/r.damflood/SWE.c
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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 <stdio.h>
#include <string.h>
#include <stdlib.h>
Expand All @@ -11,7 +28,7 @@
#include <grass/linkm.h>
#include <grass/bitmap.h>

#include "SWE.h" /* specifical dependency to the header file */
#include "SWE.h" /* Specific dependency to the header file */



Expand Down Expand Up @@ -48,9 +65,9 @@ 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 */
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is it check, check for, checking, etc that is meant for the last line?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also, there is no fall.c here

float Q, vol_res,fall, volume;
/***************************************************/
int test;
Expand All @@ -64,17 +81,17 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float



// DESCRIPTION OF METHOD (italian --> TRASLATE)
// 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

// DESCRIPTION OF METHOD
// First cycle: Calculation of new water heights at time t + 1:
echoix marked this conversation as resolved.
Show resolved Hide resolved
// - Downstream of the dam: Apply continuity equation to shallow water (?)
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Unsure of this sentence, need to remove (?) when someone understands the original meaning

// 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
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This didn't make much sense, but it was the closest that kept the same meanings of the Italian comment. Someone has a better explication?

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think swe applies where speed in plane directions are of a magnitude higher than in vertical direction, and in a lake it is questionable...

// 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++) {
Expand Down Expand Up @@ -132,8 +149,8 @@ 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
// 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];
} else if (m_v1[row][col]<0 && m_v1[row-1][col]<0) {
Expand All @@ -150,7 +167,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) {
Expand Down Expand Up @@ -180,11 +197,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);
Expand All @@ -206,8 +223,8 @@ 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("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]);
Expand All @@ -225,12 +242,12 @@ 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)
echoix marked this conversation as resolved.
Show resolved Hide resolved


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
echoix marked this conversation as resolved.
Show resolved Hide resolved
/* 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])){
Expand All @@ -255,11 +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)
echoix marked this conversation as resolved.
Show resolved Hide resolved
//*****************************************************************************
if (method==1 || method==2){

/* calcolo l'abbassamento sul lago*/
/* Calculation of the lowering of the lake*/
echoix marked this conversation as resolved.
Show resolved Hide resolved
if (num_cell!=0) {
fall = (Q * timestep-vol_res) / (num_cell * res_ew * res_ns);
} else {
Expand Down Expand Up @@ -303,14 +320,14 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float



// DESCRIPTION OF METHOD (italian --> TRASLATE)
//**********************************************************************************
// 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)
// 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)
//
// NOTA:
// u(i,j) e v (i,j) sono le velocita' medie della cella i,j
// NOTE:
// u(i,j) and v(i,j) are the average velocities of cells i,j
/*******************************************************************/
for (row = 1; row < nrows-1; row++)
{
Expand All @@ -319,7 +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 */
// 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];
Expand Down Expand Up @@ -355,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;
}
Expand All @@ -367,7 +384,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) {
Expand All @@ -385,7 +402,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) {
Expand Down Expand Up @@ -417,7 +434,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;
Expand Down Expand Up @@ -485,15 +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
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
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
/*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)
Expand All @@ -502,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");
echoix marked this conversation as resolved.
Show resolved Hide resolved
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);
Expand All @@ -525,7 +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 */
// 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];
Expand Down Expand Up @@ -576,7 +593,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) {
Expand All @@ -593,7 +610,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) {
Expand All @@ -612,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
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;
}
Expand All @@ -624,7 +641,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])) {
Expand Down Expand Up @@ -689,24 +706,24 @@ 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]); // 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]); // 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
/*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]))))>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);
Expand All @@ -726,7 +743,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]);
Expand All @@ -735,10 +752,10 @@ void shallow_water(double **m_h1,double **m_u1, double **m_v1, float **m_z,float


} else {
// tolgo h<hmin quando si svuota
// Remove h<hmin when empty
m_u2[row][col] = 0.0;
m_v2[row][col] = 0.0;
} // ciclo if (h>hmin)
} // Loop if h>hmin
}
}

Expand Down
8 changes: 4 additions & 4 deletions src/raster/r.damflood/SWE.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,10 +1,10 @@

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
/* 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
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm sure that there is something better than "In the general case", if the Italian subtleties are understood.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it is correct, but my English is not perfect ;-)

and method=3

returns void
*/
Expand Down
Loading