 Timestamp:
 Aug 28, 2007, 5:09:57 PM (16 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

anuga_core/source/anuga/shallow_water/shallow_water_ext.c
r4685 r4687 203 203 double s_min, s_max, soundspeed_left, soundspeed_right; 204 204 double denom, z; 205 double q_left_ copy[3], q_right_copy[3];205 double q_left_rotated[3], q_right_rotated[3]; 206 206 double flux_right[3], flux_left[3]; 207 207 … … 210 210 //epsilon! 211 211 212 // Copy conserved quantities to protect from modification212 // Copy conserved quantities to protect from modification 213 213 for (i=0; i<3; i++) { 214 q_left_ copy[i] = q_left[i];215 q_right_ copy[i] = q_right[i];216 } 217 218 // Align x and ymomentum with xaxis219 _rotate(q_left_ copy, n1, n2);220 _rotate(q_right_ copy, n1, n2);221 222 z = (z_left+z_right)/2; // Take average of fieldvalues223 224 // Compute speeds in xdirection225 w_left = q_left_ copy[0];214 q_left_rotated[i] = q_left[i]; 215 q_right_rotated[i] = q_right[i]; 216 } 217 218 // Align x and ymomentum with xaxis 219 _rotate(q_left_rotated, n1, n2); 220 _rotate(q_right_rotated, n1, n2); 221 222 z = (z_left+z_right)/2; // Average elevation values 223 224 // Compute speeds in xdirection 225 w_left = q_left_rotated[0]; 226 226 h_left = w_leftz; 227 uh_left = q_left_ copy[1];227 uh_left = q_left_rotated[1]; 228 228 u_left = _compute_speed(&uh_left, &h_left, epsilon, h0); 229 229 230 w_right = q_right_ copy[0];230 w_right = q_right_rotated[0]; 231 231 h_right = w_rightz; 232 uh_right = q_right_ copy[1];232 uh_right = q_right_rotated[1]; 233 233 u_right = _compute_speed(&uh_right, &h_right, epsilon, h0); 234 234 235 // Momentum in ydirection236 vh_left = q_left_ copy[2];237 vh_right = q_right_ copy[2];235 // Momentum in ydirection 236 vh_left = q_left_rotated[2]; 237 vh_right = q_right_rotated[2]; 238 238 239 239 // Limit ymomentum if necessary … … 241 241 v_right = _compute_speed(&vh_right, &h_right, epsilon, h0); 242 242 243 // Maximal and minimal wave speeds243 // Maximal and minimal wave speeds 244 244 soundspeed_left = sqrt(g*h_left); 245 245 soundspeed_right = sqrt(g*h_right); … … 252 252 253 253 254 // Flux formulas254 // Flux formulas 255 255 flux_left[0] = u_left*h_left; 256 256 flux_left[1] = u_left*uh_left + 0.5*g*h_left*h_left; … … 263 263 264 264 265 // Flux computation265 // Flux computation 266 266 denom = s_maxs_min; 267 267 if (denom == 0.0) { … … 271 271 for (i=0; i<3; i++) { 272 272 edgeflux[i] = s_max*flux_left[i]  s_min*flux_right[i]; 273 edgeflux[i] += s_max*s_min*(q_right_ copy[i]q_left_copy[i]);273 edgeflux[i] += s_max*s_min*(q_right_rotated[i]q_left_rotated[i]); 274 274 edgeflux[i] /= denom; 275 275 } 276 276 277 // Maximal wavespeed277 // Maximal wavespeed 278 278 *max_speed = max(fabs(s_max), fabs(s_min)); 279 279 280 // Rotate back280 // Rotate back 281 281 _rotate(edgeflux, n1, n2); 282 282 } … … 321 321 double s_min, s_max, soundspeed_left, soundspeed_right; 322 322 double z; 323 double q_left_ copy[3], q_right_copy[3];323 double q_left_rotated[3], q_right_rotated[3]; 324 324 325 325 double h0 = H0*H0; //This ensures a good balance when h approaches H0. … … 327 327 //Copy conserved quantities to protect from modification 328 328 for (i=0; i<3; i++) { 329 q_left_ copy[i] = q_left[i];330 q_right_ copy[i] = q_right[i];329 q_left_rotated[i] = q_left[i]; 330 q_right_rotated[i] = q_right[i]; 331 331 } 332 332 333 333 //Align x and ymomentum with xaxis 334 _rotate(q_left_ copy, n1, n2);335 _rotate(q_right_ copy, n1, n2);334 _rotate(q_left_rotated, n1, n2); 335 _rotate(q_right_rotated, n1, n2); 336 336 337 337 z = (z_left+z_right)/2; //Take average of field values 338 338 339 339 //Compute speeds in xdirection 340 w_left = q_left_ copy[0];340 w_left = q_left_rotated[0]; 341 341 h_left = w_leftz; 342 uh_left = q_left_ copy[1];342 uh_left = q_left_rotated[1]; 343 343 u_left =_compute_speed(&uh_left, &h_left, epsilon, h0); 344 344 345 w_right = q_right_ copy[0];345 w_right = q_right_rotated[0]; 346 346 h_right = w_rightz; 347 uh_right = q_right_ copy[1];347 uh_right = q_right_rotated[1]; 348 348 u_right =_compute_speed(&uh_right, &h_right, epsilon, h0); 349 349 350 350 351 351 //Momentum in ydirection 352 vh_left = q_left_ copy[2];353 vh_right = q_right_ copy[2];352 vh_left = q_left_rotated[2]; 353 vh_right = q_right_rotated[2]; 354 354 355 355 … … 702 702 703 703 PyArrayObject *h, *v, *x, *xmom, *ymom; 704 int k, i,N, k3, k6;704 int k, N, k3, k6; 705 705 double g, avg_h, zx, zy; 706 706 double x0, y0, x1, y1, x2, y2, z0, z1, z2; 707 //double epsilon; 707 708 708 709 if (!PyArg_ParseTuple(args, "dOOOOO", 709 710 &g, &h, &v, &x, 710 711 &xmom, &ymom)) { 712 //&epsilon)) { 711 713 PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: gravity could not parse input arguments"); 712 714 return NULL; … … 716 718 for (k=0; k<N; k++) { 717 719 k3 = 3*k; // base index 720 721 // Get bathymetry 722 z0 = ((double*) v > data)[k3 + 0]; 723 z1 = ((double*) v > data)[k3 + 1]; 724 z2 = ((double*) v > data)[k3 + 2]; 725 726 // Optimise for flat bed 727 // Note (Ole): This didn't produce measurable speed up. 728 // Revisit later 729 //if (fabs(z0z1)<epsilon && fabs(z1z2)<epsilon) { 730 // continue; 731 //} 732 733 // Get average depth from centroid values 734 avg_h = ((double *) h > data)[k]; 735 736 // Compute bed slope 718 737 k6 = 6*k; // base index 719 720 avg_h = 0.0; 721 for (i=0; i<3; i++) { 722 avg_h += ((double *) h > data)[k3+i]; 723 } 724 avg_h /= 3; 725 726 727 //Compute bed slope 738 728 739 x0 = ((double*) x > data)[k6 + 0]; 729 740 y0 = ((double*) x > data)[k6 + 1]; … … 734 745 735 746 736 z0 = ((double*) v > data)[k3 + 0];737 z1 = ((double*) v > data)[k3 + 1];738 z2 = ((double*) v > data)[k3 + 2];739 740 747 _gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2, &zx, &zy); 741 748 742 // Update momentum749 // Update momentum 743 750 ((double*) xmom > data)[k] += g*zx*avg_h; 744 751 ((double*) ymom > data)[k] += g*zy*avg_h; … … 1317 1324 double edgeflux[3]; // Work array for summing up fluxes 1318 1325 1319 int number_of_elements, k, i, m, n; //, j, computation_needed;1326 int number_of_elements, k, i, m, n; 1320 1327 1321 1328 int ki, nm=0, ki2; // Index shorthands … … 1420 1427 } 1421 1428 1422 1423 1424 1429 1425 1430 // Outward pointing normal vector (domain.normals[k, 2*i:2*i+2]) … … 1428 1433 normal[1] = ((double *) normals > data)[ki2+1]; 1429 1434 1435 1430 1436 // Edge flux computation (triangle k, edge i) 1431 1437 flux_function_central(ql, qr, zl, zr,
Note: See TracChangeset
for help on using the changeset viewer.