Changeset 5587 for anuga_work/development/anuga_1d/comp_flux_ext.c
 Timestamp:
 Jul 30, 2008, 5:03:47 PM (15 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

anuga_work/development/anuga_1d/comp_flux_ext.c
r5564 r5587 21 21 22 22 23 24 25 // Function to obtain speed from momentum and depth. 26 // This is used by flux functions 27 // Input parameters uh and h may be modified by this function. 28 double _compute_speed(double *uh, 29 double *h, 30 double epsilon, 31 double h0) { 32 33 double u; 34 35 if (*h < epsilon) { 36 *h = 0.0; //Could have been negative 37 u = 0.0; 38 } else { 39 u = *uh/(*h + h0/ *h); 40 } 41 42 43 // Adjust momentum to be consistent with speed 44 *uh = u * *h; 45 46 return u; 47 } 48 49 50 23 51 //Innermost flux function (using w=z+h) 24 52 int _flux_function(double *q_left, double *q_right, 25 53 double z_left, double z_right, 26 double normals, double g, double epsilon,54 double normals, double g, double epsilon, double h0, 27 55 double *edgeflux, double *max_speed) { 28 56 … … 33 61 double s_max, s_min, denom; 34 62 35 63 //printf("h0 = %f \n",h0); 36 64 ql[0] = q_left[0]; 37 65 ql[1] = q_left[1]; … … 54 82 h_left = w_leftz; 55 83 uh_left = ql[1]; 56 if (h_left < epsilon) { 57 u_left = 0.0; h_left = 0.0; 58 } 59 else { 60 u_left = uh_left/h_left; 61 } 62 84 85 u_left = _compute_speed(&uh_left, &h_left, epsilon, h0); 86 63 87 w_right = qr[0]; 64 88 h_right = w_rightz; 65 89 uh_right = qr[1]; 66 if (h_right < epsilon) { 67 u_right = 0.0; h_right = 0.0; 68 } 69 else { 70 u_right = uh_right/h_right; 71 } 90 91 u_right = _compute_speed(&uh_right, &h_right, epsilon, h0); 72 92 73 93 soundspeed_left = sqrt(g*h_left); … … 113 133 // Computational function for flux computation 114 134 double _compute_fluxes_ext(double timestep, 115 double epsilon, 116 double g, 117 long* neighbours, 118 long* neighbour_vertices, 119 double* normals, 120 double* areas, 121 double* stage_edge_values, 122 double* xmom_edge_values, 123 double* bed_edge_values, 124 double* stage_boundary_values, 125 double* xmom_boundary_values, 126 double* stage_explicit_update, 127 double* xmom_explicit_update, 128 int number_of_elements, 129 double* max_speed_array) { 130 131 double flux[2], ql[2], qr[2], edgeflux[2]; 135 double epsilon, 136 double g, 137 double h0, 138 long* neighbours, 139 long* neighbour_vertices, 140 double* normals, 141 double* areas, 142 double* stage_edge_values, 143 double* xmom_edge_values, 144 double* bed_edge_values, 145 double* stage_boundary_values, 146 double* xmom_boundary_values, 147 double* stage_explicit_update, 148 double* xmom_explicit_update, 149 int number_of_elements, 150 double* max_speed_array) { 151 152 double flux[2], ql[2], qr[2], edgeflux[2]; 132 153 double zl, zr, max_speed, normal; 133 154 int k, i, ki, n, m, nm=0; … … 159 180 160 181 normal = normals[ki]; 161 _flux_function(ql, qr, zl, zr, normal, g, epsilon, edgeflux, &max_speed);182 _flux_function(ql, qr, zl, zr, normal, g, epsilon, h0, edgeflux, &max_speed); 162 183 flux[0] = edgeflux[0]; 163 184 flux[1] = edgeflux[1]; … … 209 230 *max_speed_array; 210 231 211 double timestep, epsilon, g ;232 double timestep, epsilon, g, h0; 212 233 int number_of_elements; 213 234 214 235 // Convert Python arguments to C 215 if (!PyArg_ParseTuple(args, "ddd OOOOOOOOOOOiO",236 if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOiO", 216 237 ×tep, 217 238 &epsilon, 218 239 &g, 240 &h0, 219 241 &neighbours, 220 242 &neighbour_vertices, … … 240 262 epsilon, 241 263 g, 264 h0, 242 265 (long*) neighbours > data, 243 266 (long*) neighbour_vertices > data, … … 281 304 *max_speed_array; 282 305 283 double timestep, epsilon, g ;306 double timestep, epsilon, g, h0; 284 307 int number_of_elements; 285 308 … … 299 322 epsilon = get_python_double(domain,"epsilon"); 300 323 g = get_python_double(domain,"g"); 324 h0 = get_python_double(domain,"h0"); 301 325 302 326 … … 327 351 // the explicit update arrays 328 352 timestep = _compute_fluxes_ext(timestep, 329 epsilon, 330 g, 331 (long*) neighbours > data, 332 (long*) neighbour_vertices > data, 333 (double*) normals > data, 334 (double*) areas > data, 335 (double*) stage_vertex_values > data, 336 (double*) xmom_vertex_values > data, 337 (double*) bed_vertex_values > data, 338 (double*) stage_boundary_values > data, 339 (double*) xmom_boundary_values > data, 340 (double*) stage_explicit_update > data, 341 (double*) xmom_explicit_update > data, 342 number_of_elements, 343 (double*) max_speed_array > data); 353 epsilon, 354 g, 355 h0, 356 (long*) neighbours > data, 357 (long*) neighbour_vertices > data, 358 (double*) normals > data, 359 (double*) areas > data, 360 (double*) stage_vertex_values > data, 361 (double*) xmom_vertex_values > data, 362 (double*) bed_vertex_values > data, 363 (double*) stage_boundary_values > data, 364 (double*) xmom_boundary_values > data, 365 (double*) stage_explicit_update > data, 366 (double*) xmom_explicit_update > data, 367 number_of_elements, 368 (double*) max_speed_array > data); 344 369 345 370
Note: See TracChangeset
for help on using the changeset viewer.