- Timestamp:
- Sep 12, 2007, 1:05:19 PM (17 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/shallow_water/shallow_water_ext.c
r4730 r4731 280 280 // Flux computation 281 281 denom = s_max-s_min; 282 if (denom == 0.0) { // FIXME (Ole): Try using epsilonhere282 if (denom < epsilon) { // FIXME (Ole): Try using H0 here 283 283 for (i=0; i<3; i++) edgeflux[i] = 0.0; 284 284 *max_speed = 0.0; … … 623 623 double u, v, reduced_speed; 624 624 625 // Protect against initesimal and negative heights 626 for (k=0; k<N; k++) { 627 hc = wc[k] - zc[k]; 628 629 if (hc < minimum_allowed_height) { 625 626 // Protect against initesimal and negative heights 627 if (maximum_allowed_speed < epsilon) { 628 for (k=0; k<N; k++) { 629 hc = wc[k] - zc[k]; 630 631 if (hc < minimum_allowed_height) { 630 632 631 //Old code: Set momentum to zero and ensure h is non negative 632 //xmomc[k] = 0.0; 633 //ymomc[k] = 0.0; 634 //if (hc <= 0.0) wc[k] = zc[k]; 635 636 637 //New code: Adjust momentum to guarantee speeds are physical 638 // ensure h is non negative 639 //FIXME (Ole): This is only implemented in this C extension and 640 // has no Python equivalent 641 642 if (hc <= 0.0) { 643 wc[k] = zc[k]; 633 // Set momentum to zero and ensure h is non negative 644 634 xmomc[k] = 0.0; 645 635 ymomc[k] = 0.0; 646 } else { 647 //Reduce excessive speeds derived from division by small hc 648 //FIXME (Ole): This may be unnecessary with new slope limiters 649 //in effect. 650 651 u = xmomc[k]/hc; 652 if (fabs(u) > maximum_allowed_speed) { 653 reduced_speed = maximum_allowed_speed * u/fabs(u); 654 //printf("Speed (u) has been reduced from %.3f to %.3f\n", 655 // u, reduced_speed); 656 xmomc[k] = reduced_speed * hc; 657 } 658 659 v = ymomc[k]/hc; 660 if (fabs(v) > maximum_allowed_speed) { 661 reduced_speed = maximum_allowed_speed * v/fabs(v); 662 //printf("Speed (v) has been reduced from %.3f to %.3f\n", 663 // v, reduced_speed); 664 ymomc[k] = reduced_speed * hc; 665 } 636 if (hc <= 0.0) wc[k] = zc[k]; 637 } 638 } 639 } else { 640 641 // Protect against initesimal and negative heights 642 for (k=0; k<N; k++) { 643 hc = wc[k] - zc[k]; 644 645 if (hc < minimum_allowed_height) { 646 647 //New code: Adjust momentum to guarantee speeds are physical 648 // ensure h is non negative 649 //FIXME (Ole): This is only implemented in this C extension and 650 // has no Python equivalent 651 652 if (hc <= 0.0) { 653 wc[k] = zc[k]; 654 xmomc[k] = 0.0; 655 ymomc[k] = 0.0; 656 } else { 657 //Reduce excessive speeds derived from division by small hc 658 //FIXME (Ole): This may be unnecessary with new slope limiters 659 //in effect. 660 661 u = xmomc[k]/hc; 662 if (fabs(u) > maximum_allowed_speed) { 663 reduced_speed = maximum_allowed_speed * u/fabs(u); 664 //printf("Speed (u) has been reduced from %.3f to %.3f\n", 665 // u, reduced_speed); 666 xmomc[k] = reduced_speed * hc; 667 } 668 669 v = ymomc[k]/hc; 670 if (fabs(v) > maximum_allowed_speed) { 671 reduced_speed = maximum_allowed_speed * v/fabs(v); 672 //printf("Speed (v) has been reduced from %.3f to %.3f\n", 673 // v, reduced_speed); 674 ymomc[k] = reduced_speed * hc; 675 } 676 } 666 677 } 667 678 } … … 669 680 return 0; 670 681 } 671 672 682 673 683 … … 680 690 double cw) { 681 691 682 // Assign windfield values to momentum updates692 // Assign windfield values to momentum updates 683 693 684 694 int k;
Note: See TracChangeset
for help on using the changeset viewer.