Changeset 7154 for anuga_core/source/anuga/shallow_water
- Timestamp:
- Jun 3, 2009, 1:16:10 PM (16 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/shallow_water/shallow_water_ext.c
r7143 r7154 68 68 if (dq1>=dq2){ 69 69 if (dq1>=0.0) 70 70 *qmax=dq0+dq1; 71 71 else 72 73 72 *qmax=dq0; 73 74 74 *qmin=dq0+dq2; 75 75 if (*qmin>=0.0) *qmin = 0.0; … … 77 77 else{// dq1<dq2 78 78 if (dq2>0) 79 79 *qmax=dq0+dq2; 80 80 else 81 82 81 *qmax=dq0; 82 83 83 *qmin=dq0+dq1; 84 84 if (*qmin>=0.0) *qmin=0.0; … … 88 88 if (dq1<=dq2){ 89 89 if (dq1<0.0) 90 90 *qmin=dq0+dq1; 91 91 else 92 93 92 *qmin=dq0; 93 94 94 *qmax=dq0+dq2; 95 95 if (*qmax<=0.0) *qmax=0.0; … … 97 97 else{// dq1>dq2 98 98 if (dq2<0.0) 99 99 *qmin=dq0+dq2; 100 100 else 101 102 101 *qmin=dq0; 102 103 103 *qmax=dq0+dq1; 104 104 if (*qmax<=0.0) *qmax=0.0; … … 683 683 // FIXME: Try with this one precomputed 684 684 for (i=0; i<3; i++) { 685 685 dz = max(dz, fabs(zv[k3+i]-zc[k])); 686 686 } 687 687 } … … 711 711 712 712 if (dz > 0.0) { 713 713 alpha = max( min( alpha_balance*hmin/dz, 1.0), 0.0 ); 714 714 } else { 715 715 alpha = 1.0; // Flat bed 716 716 } 717 717 //printf("Using old style limiter\n"); … … 726 726 727 727 if (hmin < H0) { 728 alpha = 1.0; 729 for (i=0; i<3; i++) { 730 731 h_diff = hc_k - hv[i]; 732 if (h_diff <= 0) { 733 // Deep water triangle is further away from bed than 734 // shallow water (hbar < h). Any alpha will do 735 736 } else { 737 // Denominator is positive which means that we need some of the 738 // h-limited stage. 739 740 alpha = min(alpha, (hc_k - H0)/h_diff); 728 alpha = 1.0; 729 for (i=0; i<3; i++) { 730 731 h_diff = hc_k - hv[i]; 732 if (h_diff <= 0) { 733 // Deep water triangle is further away from bed than 734 // shallow water (hbar < h). Any alpha will do 735 736 } else { 737 // Denominator is positive which means that we need some of the 738 // h-limited stage. 739 740 alpha = min(alpha, (hc_k - H0)/h_diff); 741 } 742 } 743 744 // Ensure alpha in [0,1] 745 if (alpha>1.0) alpha=1.0; 746 if (alpha<0.0) alpha=0.0; 747 748 } else { 749 // Use w-limited stage exclusively in deeper water. 750 alpha = 1.0; 741 751 } 742 752 } 743 744 // Ensure alpha in [0,1] 745 if (alpha>1.0) alpha=1.0; 746 if (alpha<0.0) alpha=0.0; 747 748 } else { 749 // Use w-limited stage exclusively in deeper water. 750 alpha = 1.0; 751 } 752 } 753 754 753 754 755 755 // Let 756 756 // … … 770 770 // Momentum is balanced between constant and limited 771 771 772 772 773 773 if (alpha < 1) { 774 774 for (i=0; i<3; i++) { 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 775 776 wv[k3+i] = zv[k3+i] + (1-alpha)*hc_k + alpha*hv[i]; 777 778 // Update momentum at vertices 779 if (use_centroid_velocities == 1) { 780 // This is a simple, efficient and robust option 781 // It uses first order approximation of velocities, but retains 782 // the order used by stage. 783 784 // Speeds at centroids 785 if (hc_k > epsilon) { 786 uc = xmomc[k]/hc_k; 787 vc = ymomc[k]/hc_k; 788 } else { 789 uc = 0.0; 790 vc = 0.0; 791 } 792 793 // Vertex momenta guaranteed to be consistent with depth guaranteeing 794 // controlled speed 795 hv[i] = wv[k3+i] - zv[k3+i]; // Recompute (balanced) vertex depth 796 xmomv[k3+i] = uc*hv[i]; 797 ymomv[k3+i] = vc*hv[i]; 798 799 } else { 800 // Update momentum as a linear combination of 801 // xmomc and ymomc (shallow) and momentum 802 // from extrapolator xmomv and ymomv (deep). 803 // This assumes that values from xmomv and ymomv have 804 // been established e.g. by the gradient limiter. 805 806 // FIXME (Ole): I think this should be used with vertex momenta 807 // computed above using centroid_velocities instead of xmomc 808 // and ymomc as they'll be more representative first order 809 // values. 810 811 xmomv[k3+i] = (1-alpha)*xmomc[k] + alpha*xmomv[k3+i]; 812 ymomv[k3+i] = (1-alpha)*ymomc[k] + alpha*ymomv[k3+i]; 813 814 } 815 815 } 816 816 } … … 843 843 if (hc < minimum_allowed_height) { 844 844 845 846 847 848 845 // Set momentum to zero and ensure h is non negative 846 xmomc[k] = 0.0; 847 ymomc[k] = 0.0; 848 if (hc <= 0.0) wc[k] = zc[k]; 849 849 } 850 850 } … … 854 854 for (k=0; k<N; k++) { 855 855 hc = wc[k] - zc[k]; 856 856 857 857 if (hc < minimum_allowed_height) { 858 858 … … 866 866 } else { 867 867 //Reduce excessive speeds derived from division by small hc 868 869 868 //FIXME (Ole): This may be unnecessary with new slope limiters 869 //in effect. 870 870 871 871 u = xmomc[k]/hc; 872 873 874 875 876 877 878 872 if (fabs(u) > maximum_allowed_speed) { 873 reduced_speed = maximum_allowed_speed * u/fabs(u); 874 //printf("Speed (u) has been reduced from %.3f to %.3f\n", 875 // u, reduced_speed); 876 xmomc[k] = reduced_speed * hc; 877 } 878 879 879 v = ymomc[k]/hc; 880 881 882 883 884 885 880 if (fabs(v) > maximum_allowed_speed) { 881 reduced_speed = maximum_allowed_speed * v/fabs(v); 882 //printf("Speed (v) has been reduced from %.3f to %.3f\n", 883 // v, reduced_speed); 884 ymomc[k] = reduced_speed * hc; 885 } 886 886 } 887 887 }
Note: See TracChangeset
for help on using the changeset viewer.