- Timestamp:
- Aug 22, 2007, 4:39:07 PM (17 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/shallow_water/shallow_water_ext.c
r4631 r4677 126 126 127 127 128 void adjust_froude_number(double *uh, 129 double h, 130 double g) { 131 132 // Adjust momentum if Froude number is excessive 133 double max_froude_number = 20.0; 134 double froude_number; 135 136 //Compute Froude number (stability diagnostics) 137 froude_number = *uh/sqrt(g*h)/h; 138 139 if (froude_number > max_froude_number) { 140 printf("---------------------------------------------\n"); 141 printf("froude_number=%f (uh=%f, h=%f)\n", froude_number, *uh, h); 142 143 *uh = *uh/fabs(*uh) * max_froude_number * sqrt(g*h)*h; 144 145 froude_number = *uh/sqrt(g*h)/h; 146 printf("Adjusted froude_number=%f (uh=%f, h=%f)\n", froude_number, *uh, h); 147 printf("---------------------------------------------\n"); 148 } 149 } 150 151 152 128 153 // Function to obtain speed from momentum and depth. 129 154 // This is used by flux functions … … 135 160 136 161 double u; 162 163 //adjust_froude_number(uh, *h, 9.81); // Highly experimental and 164 // probably unneccessary 137 165 138 166 if (*h < epsilon) { … … 140 168 u = 0.0; 141 169 } else { 142 u = *uh/(*h + h0/ *h); 143 } 144 145 // printf("u=%f, h=%f\n", u, h); 170 u = *uh/(*h + h0/ *h); 171 } 172 173 174 // Adjust momentum to be consistent with speed 175 *uh = u * *h; 176 146 177 return u; 147 178 } … … 169 200 double w_left, h_left, uh_left, vh_left, u_left; 170 201 double w_right, h_right, uh_right, vh_right, u_right; 202 double v_left, v_right; 171 203 double s_min, s_max, soundspeed_left, soundspeed_right; 172 204 double denom, z; … … 175 207 176 208 double h0 = H0*H0; //This ensures a good balance when h approaches H0. 209 //But evidence suggests that h0 can be as little as 210 //epsilon! 177 211 178 212 //Copy conserved quantities to protect from modification … … 203 237 vh_right = q_right_copy[2]; 204 238 239 // Limit momentum if necessary 240 v_left =_compute_speed(&vh_left, &h_left, epsilon, h0); 241 v_right =_compute_speed(&vh_right, &h_right, epsilon, h0); 205 242 206 243 //Maximal and minimal wave speeds … … 208 245 soundspeed_right = sqrt(g*h_right); 209 246 247 210 248 s_max = max(u_left+soundspeed_left, u_right+soundspeed_right); 211 249 if (s_max < 0.0) s_max = 0.0;
Note: See TracChangeset
for help on using the changeset viewer.