Changeset 4251
 Timestamp:
 Feb 9, 2007, 4:23:15 PM (17 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

anuga_core/source/anuga/shallow_water/shallow_water_ext.c
r4239 r4251 125 125 } 126 126 127 128 // Function to obtain speed from momentum and depth. 129 // This is used by flux functions 130 // Input parameters uh and h may be modified by this function. 131 double _compute_speed(double *uh, 132 double *h, 133 double epsilon, 134 double h0) { 135 136 double u; 137 138 if (*h < epsilon) { 139 *h = 0.0; //Could have been negative 140 u = 0.0; 141 } else { 142 u = *uh/(*h + h0/ *h); 143 } 144 145 return u; 146 } 147 127 148 // Computational function for flux computation (using stage w=z+h) 128 149 int flux_function_central(double *q_left, double *q_right, … … 167 188 168 189 //Compute speeds in xdirection 169 w_left = q_left_copy[0]; // h+z190 w_left = q_left_copy[0]; 170 191 h_left = w_leftz; 171 192 uh_left = q_left_copy[1]; 172 173 if (h_left < epsilon) { 174 h_left = 0.0; //Could have been negative 175 u_left = 0.0; 176 } else { 177 u_left = uh_left/(h_left + h0/h_left); 178 } 193 u_left =_compute_speed(&uh_left, &h_left, epsilon, h0); 179 194 180 195 w_right = q_right_copy[0]; 181 196 h_right = w_rightz; 182 197 uh_right = q_right_copy[1]; 183 184 if (h_right < epsilon) { 185 h_right = 0.0; //Could have been negative 186 u_right = 0.0; 187 } else { 188 u_right = uh_right/(h_right + h0/h_right); 189 } 198 u_right =_compute_speed(&uh_right, &h_right, epsilon, h0); 190 199 191 200 //Momentum in ydirection … … 288 297 289 298 //Compute speeds in xdirection 290 w_left = q_left_copy[0]; // h+z299 w_left = q_left_copy[0]; 291 300 h_left = w_leftz; 292 301 uh_left = q_left_copy[1]; 293 294 if (h_left < epsilon) { 295 h_left = 0.0; //Could have been negative 296 u_left = 0.0; 297 } else { 298 u_left = uh_left/(h_left + h0/h_left); 299 } 302 u_left =_compute_speed(&uh_left, &h_left, epsilon, h0); 300 303 301 304 w_right = q_right_copy[0]; 302 305 h_right = w_rightz; 303 306 uh_right = q_right_copy[1]; 304 305 if (h_right < epsilon) { 306 h_right = 0.0; //Could have been negative 307 u_right = 0.0; 308 } else { 309 u_right = uh_right/(h_right + h0/h_right); 310 } 307 u_right =_compute_speed(&uh_right, &h_right, epsilon, h0); 308 311 309 312 310 //Momentum in ydirection
Note: See TracChangeset
for help on using the changeset viewer.