@@ -43,7 +43,7 @@ float velocita_breccia(int i, double h)
43
43
{
44
44
// double h;
45
45
// int i;
46
- // float g= 9.81;
46
+ // float g = 9.81;
47
47
float v ;
48
48
49
49
if (i == 1 ) {
@@ -62,7 +62,7 @@ void shallow_water(double **m_h1, double **m_u1, double **m_v1, float **m_z,
62
62
int method , int num_cell , int num_break , double t )
63
63
{
64
64
65
- /*FUNCTION VARIABLES*/
65
+ /* FUNCTION VARIABLES */
66
66
double h_dx , h_sx , h_up , h_dw , Fup , Fdw , Fdx , Fsx , Gup , Gdw , Gdx , Gsx ;
67
67
double u_sx , u_dx , v_dx , v_sx , v_up , v_dw , u_up , u_dw ;
68
68
@@ -82,23 +82,24 @@ void shallow_water(double **m_h1, double **m_u1, double **m_v1, float **m_z,
82
82
83
83
// DESCRIPTION OF METHOD
84
84
// First cycle: Calculation of new water heights at time t + 1:
85
- // - Downstream of the dam: Apply continuity equation to shallow water
86
- // In practice, the new height is evaluated through a balance
87
- // of the incoming and outgoing flows in the two main directions
88
- // - Upstream of the dam:
89
- // - In methods 1 and 2:
90
- // - The continuity equation is applied to the volume of the lake
91
- // Physically this leads to a less realistic but avoids
92
- // the oscillations that causes numerical instability.
93
- // - In the more general case the equations are applied to the whole lake
85
+ // - Downstream of the dam: Apply continuity equation to shallow water.
86
+ // In practice, the new height is evaluated through a balance
87
+ // of the incoming and outgoing flows in the two main directions
88
+ // - Upstream of the dam:
89
+ // - In methods 1 and 2:
90
+ // - The continuity equation is applied to the volume of the lake.
91
+ // Physically this leads to a less realistic but avoids the
92
+ // oscillations that causes numerical instability.
93
+ // - In the more general case, the equations are applied to the whole
94
+ // lake
94
95
95
96
for (row = 1 ; row < nrows - 1 ; row ++ ) {
96
97
for (col = 1 ; col < ncols - 1 ; col ++ ) {
97
98
if (m_lake [row ][col ] == 0 && m_DAMBREAK [row ][col ] <= 0 ) {
98
99
99
- // *******************************************/
100
- / * CONTINUITY EQUATION --> h(t+1) */
101
- // *******************************************/
100
+ /*******************************************
101
+ * CONTINUITY EQUATION --> h(t+1)
102
+ *******************************************/
102
103
// x direction
103
104
// right intercell
104
105
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,
168
169
}
169
170
F = Fdx - Fsx ;
170
171
171
- // dGup =m_v1[row][col] * m_h1[row][col]; y direction
172
- // intercell up
172
+ /* dGup = m_v1[row][col] * m_h1[row][col]; */
173
+ // y direction
174
+ // intercell up
173
175
if (m_v1 [row ][col ] > 0 && m_v1 [row - 1 ][col ] > 0 ) {
174
176
Gup = m_v1 [row ][col ] * m_h1 [row ][col ];
175
177
}
@@ -242,29 +244,48 @@ void shallow_water(double **m_h1, double **m_u1, double **m_v1, float **m_z,
242
244
m_h2 [row ][col ] = m_h1 [row ][col ] - timestep / res_ew * F -
243
245
timestep / res_ns * G ;
244
246
245
- /*if ((row==20||row==21||row==22||row==23)&&(col==18||col==19)){
246
- printf("EQ. CONTINUITY --> row:%d, col:%d\n)",row,
247
- col);
248
- 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]);
249
- printf("m_h1[row][col+1]:%f,m_h1[row][col-1]:%f,m_h1[row+1][col]:%f,
250
- m_h1[row-1][col]:%f\n",m_h1[row][col+1],m_h1[row][col-1],m_h1[row+1][col],
251
- m_h1[row-1][col]); printf("h_dx:%f, h_sx:%f, h_up%f,
252
- h_dw:%f\n",h_dx, h_sx, h_up, h_dw);
253
- printf("m_u1[row][col+1]:%f,m_u1[row][col-1]:%f,m_v1[row+1][col]:%f,
254
- m_v1[row-1][col]:%f\n",m_u1[row][col+1],m_u1[row][col-1],m_v1[row+1][col],
255
- m_v1[row-1][col]); printf("v_up: %f,v_dw:%f,u_dx:%f,u_sx:%f
256
- \n",v_up, v_dw, u_dx, u_sx); printf("Fdx: %f,Fsx: %f, F: %f,
257
- Gup:%f, Gdw:%f, G: %f\n",Fdx, Fsx, F,Gup, Gdw, G);
258
- printf("m_h2(row,col): %f\n \n", m_h2[row][col]);
259
- }*/
247
+ /* if ((row == 20 || row == 21 || row == 22 || row == 23) &&
248
+ (col == 18 || col == 19)) {
249
+ printf("EQ. CONTINUITY --> row:%d, col:%d\n)", row, col);
250
+ printf("m_h1[row][col]:%f, "
251
+ "m_u1[row][col]:%f, "
252
+ "m_v1[row][col]:%f",
253
+ m_h1[row][col], m_u1[row][col], m_v1[row][col]);
254
+ printf("m_h1[row][col+1]:%f, "
255
+ "m_h1[row][col-1]:%f, "
256
+ "m_h1[row+1][col]:%f, "
257
+ "m_h1[row-1][col]:% f\n",
258
+ m_h1[row][col + 1], m_h1[row][col - 1],
259
+ m_h1[row + 1][col], m_h1[row - 1][col]);
260
+ printf("h_dx:%f, h_sx:%f, h_up%f, h_dw:%f\n", h_dx, h_sx,
261
+ h_up, h_dw);
262
+ printf("m_u1[row][col+1]:%f, "
263
+ "m_u1[row][col-1]:%f, "
264
+ "m_v1[row+1][col]:%f, "
265
+ "m_v1[row-1][col]:%f\n",
266
+ m_u1[row][col + 1], m_u1[row][col - 1],
267
+ m_v1[row + 1][col], m_v1[row - 1][col]);
268
+ printf("v_up: %f, v_dw:%f, "
269
+ "u_dx:%f, u_sx:%f\n",
270
+ v_up, v_dw, u_dx, u_sx);
271
+ printf("Fdx: %f,Fsx: %f, F: %f, "
272
+ "Gup:%f, Gdw:%f, G: %f\n",
273
+ Fdx, Fsx, F, Gup, Gdw, G);
274
+ printf("m_h2(row,col): %f\n \n", m_h2[row][col]);
275
+ } */
260
276
261
- /*if( (row==1 || row==(nrows-2) || col==1 || col==(ncols-2)) &&
262
- (m_v2[1][col]>0 || m_v2[nrows-2][col]<0 || m_u1[row][1]<0 ||
263
- m_u1[row][ncols-2]>0 )){ if (warn==0){ G_warning("At the time %.3f the
264
- computational region is smaller than inundation",t); warn=1;
265
- G_message("warn=%d",warn);
266
- }
267
- }*/
277
+ /* if (((row == 1 || row == (nrows - 2)) ||
278
+ (col == 1 || col == (ncols - 2))) &&
279
+ (m_v2[1][col] > 0 || m_v2[nrows - 2][col] < 0 ||
280
+ m_u1[row][1] < 0 || m_u1[row][ncols - 2] > 0)) {
281
+ if (warn == 0) {
282
+ G_warning("At the time %.3f the computational region "
283
+ "is smaller than inundation",
284
+ t);
285
+ warn = 1;
286
+ G_message("warn=%d", warn);
287
+ }
288
+ } */
268
289
269
290
if (m_h2 [row ][col ] < 0 ) {
270
291
/*G_warning("At the time %f h is lesser than 0
@@ -301,7 +322,8 @@ void shallow_water(double **m_h1, double **m_u1, double **m_v1, float **m_z,
301
322
302
323
if (method == 1 || method == 2 ) {
303
324
//*******************************************************************
304
- // Calculation of flow rate Q coming out of the lake only in the case of spillway Hp (hypothesis)
325
+ // Calculation of flow rate Q coming out of the lake only in the
326
+ // case of spillway Hp (hypothesis)
305
327
/* Hypothesis: method 1 or 2 */
306
328
if (m_DAMBREAK [row ][col ] > 0 ) {
307
329
if ((m_z [row ][col ] + m_h1 [row ][col ]) >
@@ -670,12 +692,12 @@ void shallow_water(double **m_h1, double **m_u1, double **m_v1, float **m_z,
670
692
(m_z [row ][col + 1 ] + m_h2 [row ][col + 1 ]))
671
693
m_u2 [row ][col ] = velocita_breccia (
672
694
method ,
673
- m_h2 [row ][col ]); //velocity on the spillway
695
+ m_h2 [row ][col ]); // velocity on the spillway
674
696
else if ((m_z [row ][col ] + m_h2 [row ][col ]) >
675
697
(m_z [row ][col - 1 ] + m_h2 [row ][col - 1 ]))
676
698
m_u2 [row ][col ] = - velocita_breccia (
677
699
method ,
678
- m_h2 [row ][col ]); //velocity on the spillway
700
+ m_h2 [row ][col ]); // velocity on the spillway
679
701
else
680
702
m_u2 [row ][col ] = 0.0 ;
681
703
}
0 commit comments