Changeset 4687
- Timestamp:
- Aug 28, 2007, 5:09:57 PM (17 years ago)
- Location:
- anuga_core/source/anuga
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/anuga_config.py
r2528 r4687 3 3 4 4 epsilon = 1.0e-12 5 6 5 default_boundary_tag = 'exterior' -
anuga_core/source/anuga/shallow_water/shallow_water_domain.py
r4685 r4687 1500 1500 ymom = domain.quantities['ymomentum'].explicit_update 1501 1501 1502 Stage = domain.quantities['stage'] 1503 Elevation = domain.quantities['elevation'] 1504 h = Stage.edge_values - Elevation.edge_values 1505 v = Elevation.vertex_values 1502 stage = domain.quantities['stage'] 1503 elevation = domain.quantities['elevation'] 1504 1505 h = stage.centroid_values - elevation.centroid_values 1506 z = elevation.vertex_values 1506 1507 1507 1508 x = domain.get_vertex_coordinates() 1508 1509 g = domain.g 1509 1510 1510 1511 1511 1512 from shallow_water_ext import gravity 1512 gravity(g, h, v, x, xmom, ymom)1513 gravity(g, h, z, x, xmom, ymom) #, 1.0e-6) 1513 1514 1514 1515 -
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 y-momentum with x-axis219 _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 x-direction225 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 y-momentum with x-axis 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 x-direction 225 w_left = q_left_rotated[0]; 226 226 h_left = w_left-z; 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_right-z; 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 y-direction236 vh_left = q_left_ copy[2];237 vh_right = q_right_ copy[2];235 // Momentum in y-direction 236 vh_left = q_left_rotated[2]; 237 vh_right = q_right_rotated[2]; 238 238 239 239 // Limit y-momentum 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_max-s_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 y-momentum with x-axis 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 x-direction 340 w_left = q_left_ copy[0];340 w_left = q_left_rotated[0]; 341 341 h_left = w_left-z; 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_right-z; 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 y-direction 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(z0-z1)<epsilon && fabs(z1-z2)<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.