Changeset 9260
Legend:
- Unmodified
- Added
- Removed
-
trunk/anuga_core/source/anuga/shallow_water/shallow_water_domain.py
r9257 r9260 250 250 self.riverwallData=anuga.structures.riverwall.RiverWall(self) 251 251 252 # Keep track of the fluxes through the boundaries253 # Only works for DE algorithms at present252 ## Keep track of the fluxes through the boundaries 253 ## Only works for DE algorithms at present 254 254 max_time_substeps=3 # Maximum number of substeps supported by any timestepping method 255 # boundary_flux_sum holds boundary fluxes on each sub-step 255 # boundary_flux_sum holds boundary fluxes on each sub-step [unused substeps = 0.] 256 256 self.boundary_flux_sum=num.array([0.]*max_time_substeps) 257 257 from anuga.operators.boundary_flux_integral_operator import boundary_flux_integral_operator … … 260 260 self.call=1 261 261 262 # Work arrays [avoid allocate statements in compute_fluxes or extrapolate_second_order] 263 self.edge_flux_work=num.zeros(len(self.edge_coordinates[:,0])*3) 264 self.pressuregrad_work=num.zeros(len(self.edge_coordinates[:,0])) 265 self.x_centroid_work=num.zeros(len(self.edge_coordinates[:,0])/3) 266 self.y_centroid_work=num.zeros(len(self.edge_coordinates[:,0])/3) 262 267 263 268 def set_defaults(self): … … 1590 1595 self.riverwallData.hydraulic_properties_rowIndex, 1591 1596 int(self.riverwallData.ncol_hydraulic_properties), 1592 self.riverwallData.hydraulic_properties) 1597 self.riverwallData.hydraulic_properties, 1598 self.edge_flux_work, 1599 self.pressuregrad_work) 1593 1600 1594 1601 self.flux_timestep = flux_timestep … … 1696 1703 height.vertex_values, 1697 1704 int(self.optimise_dry_cells), 1698 int(self.extrapolate_velocity_second_order)) 1705 int(self.extrapolate_velocity_second_order), 1706 self.x_centroid_work, 1707 self.y_centroid_work) 1699 1708 else: 1700 1709 # Code for original method … … 2411 2420 message = '---------------------------\n' 2412 2421 message += 'Volumetric balance report:\n' 2422 message += 'Note: Boundary fluxes are not exact\n' 2423 message += 'See get_boundary_flux_integral for exact computation\n' 2413 2424 message += '--------------------------\n' 2414 2425 message += 'Total boundary inflow [m^3/s]: %.2f\n' % total_boundary_inflow -
trunk/anuga_core/source/anuga/shallow_water/swDE1_domain_ext.c
r9257 r9260 448 448 double newFlux; 449 449 // Following constants control the 'blending' with the shallow water solution 450 // They are now user-defined 450 // They are now user-defined 451 451 //double s1=0.9; // At this submergence ratio, begin blending with shallow water solution 452 452 //double s2=0.95; // At this submergence ratio, completely use shallow water solution … … 474 474 } 475 475 476 //printf("%e, %e \n", rw, edgeflux[0]);477 478 476 if( (hdRat<s2) & (hdWrRat< h2) ){ 479 477 // Rescale the edge fluxes so that the mass flux = desired flux … … 489 487 // Weighted average constants to transition to shallow water eqn flow 490 488 w1=min( max(hdRat-s1,0.)/(s2-s1), 1.0); 491 //w2=1.0-w1; 492 // 489 493 490 // Adjust again when the head is too deep relative to the weir height 494 491 w2=min( max(hdWrRat-h1,0.)/(h2-h1), 1.0); 495 //w2=1.0-w1;496 492 497 493 newFlux=(rw*(1.0-w1)+w1*edgeflux[0])*(1.0-w2) + w2*edgeflux[0]; … … 505 501 edgeflux[2]*=scaleFlux; 506 502 }else{ 507 // Can't divide by edgeflux 503 // Can't divide by edgeflux, so enforce 'newFlux' directly 504 // 508 505 // FIXME: This is fluxing mass but not momentum 509 506 // It might be ok, but,......... … … 522 519 } 523 520 } 524 525 //printf("%e, %e, %e, %e , %e, %e \n", weir_height, h_left, h_right, edgeflux[0], w1, w2);526 521 527 522 return 0; … … 563 558 double* bed_centroid_values, 564 559 double* height_centroid_values, 565 //double* bed_vertex_values,566 560 long* edge_flux_type, 567 561 double* riverwall_elevation, 568 562 long* riverwall_rowIndex, 569 563 int ncol_riverwall_hydraulic_properties, 570 double* riverwall_hydraulic_properties) { 564 double* riverwall_hydraulic_properties, 565 double* edge_flux_work, 566 double* pressuregrad_work) { 571 567 // Local variables 572 568 double max_speed, length, inv_area, zl, zr; 573 //double h0 = H0*H0;//H0*H0;//H0*0.1;//H0*H0;//H0*H0; // This ensures a good balance when h approaches H0.574 569 double h_left, h_right, z_half ; // For andusse scheme 575 570 // FIXME: limiting_threshold is not used for DE1 576 double limiting_threshold = 10*H0;//10 * H0 ;//H0; //10 * H0; // Avoid applying limiter below this 577 // threshold for performance reasons. 578 // See ANUGA manual under flux limiting 571 double limiting_threshold = 10*H0; 572 // 579 573 int k, i, m, n,j, ii; 580 574 int ki, nm = 0, ki2,ki3, nm3; // Index shorthands 581 //int num_hydraulic_properties=1; // FIXME: Move to function call582 575 // Workspace (making them static actually made function slightly slower (Ole)) 583 576 double ql[3], qr[3], edgeflux[3]; // Work array for summing up fluxes … … 592 585 double speed_max_last, vol, smooth, local_speed, weir_height; 593 586 594 double *edgeflux_store, *pressuregrad_store;595 596 edgeflux_store = malloc(number_of_elements*9*sizeof(double));597 pressuregrad_store = malloc(number_of_elements*3*sizeof(double));598 599 587 call++; // Flag 'id' of flux calculation for this timestep 600 588 … … 604 592 memset((char*) xmom_explicit_update, 0, number_of_elements * sizeof (double)); 605 593 memset((char*) ymom_explicit_update, 0, number_of_elements * sizeof (double)); 594 memset((char*) edge_flux_work, 0, 9*number_of_elements * sizeof (double)); 595 memset((char*) pressuregrad_work, 0, 3*number_of_elements * sizeof (double)); 606 596 607 597 local_timestep=timestep; … … 680 670 z_half=max(riverwall_elevation[RiverWall_count-1], z_half) ; 681 671 682 //if(min(ql[0], qr[0]) < z_half){683 // // Since there is a wall blocking the flow connection, use first order extrapolation for this edge684 // ql[0]=stage_centroid_values[k];685 // ql[1]=xmom_centroid_values[k];686 // ql[2]=ymom_centroid_values[k];687 // hle=hc;688 // zl=zc;689 690 // if(n>=0){691 // qr[0]=stage_centroid_values[n];692 // qr[1]=xmom_centroid_values[n];693 // qr[2]=ymom_centroid_values[n];694 // hre=hc_n;695 // zr = zc_n;696 // }else{697 // hre=hc;698 // zr = zc;699 // }700 // // Re-set central bed to riverwall elevation701 // z_half=max(riverwall_elevation[RiverWall_count-1], max(zl, zr)) ;702 //}703 704 672 } 705 673 … … 720 688 // Force weir discharge to match weir theory 721 689 if(edge_flux_type[ki]==1){ 722 //printf("%e \n", z_half);723 690 weir_height=max(riverwall_elevation[RiverWall_count-1]-min(zl,zr), 0.); // Reference weir height 724 //weir_height=max(z_half-max(zl,zr), 0.); // Reference weir height725 691 // If the weir height is zero, avoid the weir compuattion entirely 726 692 if(weir_height>0.){ … … 728 694 ii=riverwall_rowIndex[RiverWall_count-1]*ncol_riverwall_hydraulic_properties; 729 695 Qfactor=riverwall_hydraulic_properties[ii]; 730 //printf("%e \n", Qfactor);731 696 // Get s1, submergence ratio at which we start blending with the shallow water solution 732 697 ii+=1; … … 742 707 h2=riverwall_hydraulic_properties[ii]; 743 708 744 //adjust_edgeflux_with_weir(edgeflux, h_left, h_right, g,745 // weir_height, Qfactor,746 // s1, s2, h1, h2);747 748 709 // Use first-order h's for weir -- as the 'upstream/downstream' heads are 749 710 // measured away from the weir itself … … 759 720 s1, s2, h1, h2); 760 721 // NOTE: Should perhaps also adjust the wave-speed here?? Small chance of instability?? 761 //printf("%e \n", edgeflux[0]);762 722 } 763 723 } … … 770 730 771 731 772 //// Don't allow an outward advective flux if the cell centroid stage 773 //// is < the edge value. Is this important (??) 732 //// Don't allow an outward advective flux if the cell centroid 733 //// stage is < the edge value. Is this important (??). Seems not 734 //// to be with DE algorithms 774 735 //if((hc<H0) && edgeflux[0] > 0.){ 775 736 // edgeflux[0] = 0.; … … 788 749 //} 789 750 790 edge flux_store[ki3 + 0 ] = -edgeflux[0];791 edge flux_store[ki3 + 1 ] = -edgeflux[1];792 edge flux_store[ki3 + 2 ] = -edgeflux[2];751 edge_flux_work[ki3 + 0 ] = -edgeflux[0]; 752 edge_flux_work[ki3 + 1 ] = -edgeflux[1]; 753 edge_flux_work[ki3 + 2 ] = -edgeflux[2]; 793 754 794 755 // bedslope_work contains all gravity related terms 795 756 bedslope_work=length*(-g*0.5*(h_left*h_left - hle*hle -(hle+hc)*(zl-zc))+pressure_flux); 796 757 797 pressuregrad_ store[ki]=bedslope_work;758 pressuregrad_work[ki]=bedslope_work; 798 759 799 760 already_computed_flux[ki] = call; // #k Done … … 802 763 if (n >= 0) { 803 764 804 edge flux_store[nm3 + 0 ] = edgeflux[0];805 edge flux_store[nm3 + 1 ] = edgeflux[1];806 edge flux_store[nm3 + 2 ] = edgeflux[2];765 edge_flux_work[nm3 + 0 ] = edgeflux[0]; 766 edge_flux_work[nm3 + 1 ] = edgeflux[1]; 767 edge_flux_work[nm3 + 2 ] = edgeflux[2]; 807 768 bedslope_work=length*(-g*0.5*(h_right*h_right-hre*hre-(hre+hc_n)*(zr-zc_n))+pressure_flux); 808 pressuregrad_ store[nm]=bedslope_work;769 pressuregrad_work[nm]=bedslope_work; 809 770 810 771 already_computed_flux[nm] = call; // #n Done … … 815 776 // steps, NOT within them (since a constant timestep is used within 816 777 // each rk2/rk3 sub-step) 817 if ((tri_full_flag[k] == 1) & ( (call-2)%timestep_fluxcalls==0)) {778 if ((tri_full_flag[k] == 1) & (substep_count==0)) { 818 779 819 780 speed_max_last=max(speed_max_last, max_speed); … … 833 794 } 834 795 835 // Keep track of maximal speeds836 //max_speed_array[k] = max(max_speed_array[k],max_speed);837 838 796 } // End edge i (and neighbour n) 839 797 // Keep track of maximal speeds … … 843 801 } // End triangle k 844 802 845 //// GD HACK846 803 //// Limit edgefluxes, for mass conservation near wet/dry cells 847 804 //// This doesn't seem to be needed anymore … … 855 812 // outgoing_mass_edges=0.0; 856 813 // for(useint=0; useint<3; useint++){ 857 // if(edge flux_store[3*(3*k+useint)]< 0.){814 // if(edge_flux_work[3*(3*k+useint)]< 0.){ 858 815 // //outgoing_mass_edges+=1.0; 859 // outgoing_mass_edges+=(edge flux_store[3*(3*k+useint)]);816 // outgoing_mass_edges+=(edge_flux_work[3*(3*k+useint)]); 860 817 // } 861 818 // } … … 871 828 // // total_outgoing_flux <= cell volume = Area_triangle*hc 872 829 // vol=areas[k]*hc; 873 // if((edge flux_store[ki3]< 0.0) && (-outgoing_mass_edges> vol)){830 // if((edge_flux_work[ki3]< 0.0) && (-outgoing_mass_edges> vol)){ 874 831 // 875 832 // // This bound could be improved (e.g. we could actually sum the … … 880 837 // tmp = vol/(-(outgoing_mass_edges)) ; 881 838 // if(tmp< 1.0){ 882 // edge flux_store[ki3+0]*=tmp;883 // edge flux_store[ki3+1]*=tmp;884 // edge flux_store[ki3+2]*=tmp;839 // edge_flux_work[ki3+0]*=tmp; 840 // edge_flux_work[ki3+1]*=tmp; 841 // edge_flux_work[ki3+2]*=tmp; 885 842 886 843 // // Compute neighbour edge index … … 889 846 // nm = 3*n + neighbour_edges[ki]; 890 847 // nm3 = nm*3; 891 // edge flux_store[nm3+0]*=tmp;892 // edge flux_store[nm3+1]*=tmp;893 // edge flux_store[nm3+2]*=tmp;848 // edge_flux_work[nm3+0]*=tmp; 849 // edge_flux_work[nm3+1]*=tmp; 850 // edge_flux_work[nm3+2]*=tmp; 894 851 // } 895 852 // } … … 897 854 // } 898 855 // } 899 900 ////printf("%e \n", edgeflux_store[3*30*3]);901 856 902 857 // Now add up stage, xmom, ymom explicit updates … … 916 871 // Option to limit advective fluxes 917 872 //if(hc > H0){ 918 stage_explicit_update[k] += edge flux_store[ki3+0];919 xmom_explicit_update[k] += edge flux_store[ki3+1];920 ymom_explicit_update[k] += edge flux_store[ki3+2];873 stage_explicit_update[k] += edge_flux_work[ki3+0]; 874 xmom_explicit_update[k] += edge_flux_work[ki3+1]; 875 ymom_explicit_update[k] += edge_flux_work[ki3+2]; 921 876 //}else{ 922 // stage_explicit_update[k] += edge flux_store[ki3+0];877 // stage_explicit_update[k] += edge_flux_work[ki3+0]; 923 878 //} 924 879 … … 930 885 // boundary_flux_sum is an array with length = timestep_fluxcalls 931 886 // For each sub-step, we put the boundary flux sum in. 932 boundary_flux_sum[substep_count] += edge flux_store[ki3];887 boundary_flux_sum[substep_count] += edge_flux_work[ki3]; 933 888 } 889 934 890 // GD HACK 935 891 // Compute bed slope term 936 892 //if(hc > H0){ 937 xmom_explicit_update[k] -= normals[ki2]*pressuregrad_ store[ki];938 ymom_explicit_update[k] -= normals[ki2+1]*pressuregrad_ store[ki];893 xmom_explicit_update[k] -= normals[ki2]*pressuregrad_work[ki]; 894 ymom_explicit_update[k] -= normals[ki2+1]*pressuregrad_work[ki]; 939 895 //}else{ 940 896 // xmom_explicit_update[k] *= 0.; … … 953 909 } // end cell k 954 910 955 // Hack to ensure we only update the timestep on the first call within each rk2/rk3 step911 // Ensure we only update the timestep on the first call within each rk2/rk3 step 956 912 if(substep_count==0) timestep=local_timestep; 957 913 958 free(edgeflux_store);959 free(pressuregrad_store);960 961 914 return timestep; 962 915 } … … 1045 998 1046 999 int limit_gradient(double *dqv, double qmin, double qmax, double beta_w){ 1047 // Given provisional jumps dqv from the FV triangle centroid to its 1048 // vertices and jumps qmin (qmax) between the centroid of the FV1049 // triangle and the minimum (maximum) of the values at the centroid of1050 // the FV triangle and the auxiliary triangle vertices,1051 // calculate a multiplicative factor phi by which the provisional1052 // vertex jumps are to belimited1000 // Given provisional jumps dqv from the FV triangle centroid to its 1001 // vertices/edges, and jumps qmin (qmax) between the centroid of the FV 1002 // triangle and the minimum (maximum) of the values at the auxiliary triangle 1003 // vertices (which are centroids of neighbour mesh triangles), calculate a 1004 // multiplicative factor phi by which the provisional vertex jumps are to be 1005 // limited 1053 1006 1054 1007 int i; … … 1111 1064 double* height_vertex_values, 1112 1065 int optimise_dry_cells, 1113 int extrapolate_velocity_second_order) { 1066 int extrapolate_velocity_second_order, 1067 double* x_centroid_work, 1068 double* y_centroid_work) { 1114 1069 1115 1070 // Local variables … … 1122 1077 double dk, dv0, dv1, dv2, de[3], demin, dcmax, r0scale, vel_norm, l1, l2, a_tmp, b_tmp, c_tmp,d_tmp; 1123 1078 1124 double *xmom_centroid_store, *ymom_centroid_store; // , *min_elevation_edgevalue, *max_elevation_edgevalue; 1125 //double *cell_wetness_scale; 1126 //int *count_wet_neighbours; 1127 1128 // Use malloc to avoid putting these variables on the stack, which can cause 1129 // segfaults in large model runs 1130 xmom_centroid_store = malloc(number_of_elements*sizeof(double)); 1131 ymom_centroid_store = malloc(number_of_elements*sizeof(double)); 1132 //min_elevation_edgevalue = malloc(number_of_elements*sizeof(double)); 1133 //max_elevation_edgevalue = malloc(number_of_elements*sizeof(double)); 1134 //cell_wetness_scale = malloc(number_of_elements*sizeof(double)); 1135 //count_wet_neighbours = malloc(number_of_elements*sizeof(int)); 1136 1137 // 1138 // Compute some useful statistics on wetness/dryness 1139 // 1140 //for(k=0; k<number_of_elements;k++){ 1141 // //cell_wetness_scale[k] = 0.; 1142 // //// Check if cell k is wet 1143 // ////if(stage_centroid_values[k] > elevation_centroid_values[k]){ 1144 // //if(stage_centroid_values[k] > elevation_centroid_values[k] + 1.0*minimum_allowed_height){ 1145 // ////if(stage_centroid_values[k] > max_elevation_edgevalue[k] + minimum_allowed_height+epsilon){ 1146 // // cell_wetness_scale[k] = 1.; 1147 // //} 1148 1149 // //min_elevation_edgevalue[k] = min(elevation_edge_values[3*k], 1150 // // min(elevation_edge_values[3*k+1], 1151 // // elevation_edge_values[3*k+2])); 1152 // //max_elevation_edgevalue[k] = max(elevation_edge_values[3*k], 1153 // // max(elevation_edge_values[3*k+1], 1154 // // elevation_edge_values[3*k+2])); 1155 //} 1156 1157 //// Alternative 'PROTECT' step 1158 //for(k=0; k<number_of_elements;k++){ 1159 // //if((cell_wetness_scale[k]==0. ) ){ 1160 // // xmom_centroid_values[k]=0.; 1161 // // ymom_centroid_values[k]=0.; 1162 // // //xmom_centroid_values[k]*=0.9; 1163 // // //ymom_centroid_values[k]*=0.9; 1164 1165 // //} 1166 1167 // 1168 // ////// Try some Froude-number limiting in shallow depths 1169 // //dpth=max(stage_centroid_values[k] - elevation_centroid_values[k], 0.0); 1170 // //// 1171 // //if(dpth< max(max_elevation_edgevalue[k]-min_elevation_edgevalue[k],10*minimum_allowed_height)){ 1172 // // // momnorm = momentum^2 1173 // // momnorm=(xmom_centroid_values[k]*xmom_centroid_values[k]+ 1174 // // ymom_centroid_values[k]*ymom_centroid_values[k]); 1175 // // 1176 // // // vel^2 < constant*g*dpth [-- then multiply both sides by dpth^2] 1177 // // if(momnorm > 4*9.81*dpth*dpth*dpth){ 1178 // // // Down-scale momentum so that Froude number < constant 1179 // // tmp=sqrt((4*9.81*dpth*dpth*dpth)/momnorm); 1180 // // xmom_centroid_values[k] *=tmp; 1181 // // ymom_centroid_values[k] *=tmp; 1182 // // } 1183 // //} 1184 //} 1185 1079 1080 memset((char*) x_centroid_work, 0, number_of_elements * sizeof (double)); 1081 memset((char*) y_centroid_work, 0, number_of_elements * sizeof (double)); 1082 1186 1083 1187 1084 if(extrapolate_velocity_second_order==1){ … … 1195 1092 dk = height_centroid_values[k]; 1196 1093 if(dk>minimum_allowed_height){ 1197 x mom_centroid_store[k] = xmom_centroid_values[k];1094 x_centroid_work[k] = xmom_centroid_values[k]; 1198 1095 xmom_centroid_values[k] = xmom_centroid_values[k]/dk; 1199 1096 1200 y mom_centroid_store[k] = ymom_centroid_values[k];1097 y_centroid_work[k] = ymom_centroid_values[k]; 1201 1098 ymom_centroid_values[k] = ymom_centroid_values[k]/dk; 1202 1099 }else{ 1203 x mom_centroid_store[k] = 0.;1100 x_centroid_work[k] = 0.; 1204 1101 xmom_centroid_values[k] = 0.; 1205 y mom_centroid_store[k] = 0.;1102 y_centroid_work[k] = 0.; 1206 1103 ymom_centroid_values[k] = 0.; 1207 1104 … … 1224 1121 (height_centroid_values[k1] < minimum_allowed_height | k1==k) & 1225 1122 (height_centroid_values[k2] < minimum_allowed_height | k2==k)){ 1226 x mom_centroid_store[k] = 0.;1123 x_centroid_work[k] = 0.; 1227 1124 xmom_centroid_values[k] = 0.; 1228 y mom_centroid_store[k] = 0.;1125 y_centroid_work[k] = 0.; 1229 1126 ymom_centroid_values[k] = 0.; 1230 1127 … … 1313 1210 k2 = surrogate_neighbours[k3 + 2]; 1314 1211 1315 //if(cell_wetness_scale[k0]==0 && cell_wetness_scale[k1]==0 && cell_wetness_scale[k2]==0){1316 // xmom_centroid_store[k]=0.;1317 // ymom_centroid_store[k]=0.;1318 // xmom_centroid_values[k] = 0.;1319 // ymom_centroid_values[k] = 0.;1320 //}1321 1322 // Test to see whether we accept the surrogate neighbours1323 // Note that if ki is replaced with k in more than 1 neighbour, then the1324 // triangle area will be zero, and a first order extrapolation will be1325 // used1326 // FIXME: Remove cell_wetness_scale if you don't need it1327 //if( (cell_wetness_scale[k2]==0.0 && stage_centroid_values[k]<max_elevation_edgevalue[k]+100.)){1328 // k2 = k ;1329 //}1330 //if((cell_wetness_scale[k0]==0.0 && stage_centroid_values[k]<max_elevation_edgevalue[k])+100.){1331 // k0 = k ;1332 //}1333 //if((cell_wetness_scale[k1]==0.0 && stage_centroid_values[k]<max_elevation_edgevalue[k])+100.){1334 // k1 = k ;1335 //}1336 1337 // Take note if the max neighbour bed elevation is greater than the min1338 // neighbour stage -- suggests a 'steep' bed relative to the flow1339 //bedmax = max(elevation_centroid_values[k],1340 // max(elevation_centroid_values[k0],1341 // max(elevation_centroid_values[k1],1342 // elevation_centroid_values[k2])));1343 //bedmin = min(elevation_centroid_values[k],1344 // min(elevation_centroid_values[k0],1345 // min(elevation_centroid_values[k1],1346 // elevation_centroid_values[k2])));1347 //stagemin = min(max(stage_centroid_values[k], elevation_centroid_values[k]),1348 // min(max(stage_centroid_values[k0], elevation_centroid_values[k0]),1349 // min(max(stage_centroid_values[k1], elevation_centroid_values[k1]),1350 // max(stage_centroid_values[k2], elevation_centroid_values[k2]))));1351 1352 1212 // Get the auxiliary triangle's vertex coordinates 1353 1213 // (really the centroids of neighbouring triangles) … … 1377 1237 1378 1238 //// Treat triangles with no neighbours (area2 <=0.) 1379 if ((area2 <= 0.)) //|( cell_wetness_scale[k2]==0. && cell_wetness_scale[k1]==0. && cell_wetness_scale[k0]==0.))//|| (cell_wetness_scale[k]==0.)) //|(count_wet_neighbours[k]==0))1239 if ((area2 <= 0.)) 1380 1240 { 1381 1241 … … 1400 1260 //xmom_centroid_values[k]=0.; 1401 1261 //ymom_centroid_values[k]=0.; 1402 //x mom_centroid_store[k] = 0.;1403 //y mom_centroid_store[k] = 0.;1262 //x_centroid_work[k] = 0.; 1263 //y_centroid_work[k] = 0.; 1404 1264 1405 1265 xmom_edge_values[k3] = xmom_centroid_values[k]; … … 1414 1274 1415 1275 // Calculate heights of neighbouring cells 1416 hc = height_centroid_values[k]; //stage_centroid_values[k] - elevation_centroid_values[k];1417 h0 = height_centroid_values[k0]; // - elevation_centroid_values[k0];1418 h1 = height_centroid_values[k1]; // - elevation_centroid_values[k1];1419 h2 = height_centroid_values[k2]; // - elevation_centroid_values[k2];1276 hc = height_centroid_values[k]; 1277 h0 = height_centroid_values[k0]; 1278 h1 = height_centroid_values[k1]; 1279 h2 = height_centroid_values[k2]; 1420 1280 1421 1281 hmin = min(min(h0, min(h1, h2)), hc); … … 1428 1288 // but is also more 'artefacty' in important cases (tendency for high velocities, etc). 1429 1289 // 1430 //hfactor=1.0;1431 1290 a_tmp=0.3; // Highest depth ratio with hfactor=1 1432 1291 b_tmp=0.1; // Highest depth ratio with hfactor=0 1433 1292 c_tmp=1.0/(a_tmp-b_tmp); 1434 1293 d_tmp= 1.0-(c_tmp*a_tmp); 1294 // So hfactor = depth_ratio*(c_tmp) + d_tmp, but is clipped between 0 and 1. 1435 1295 hfactor= max(0., min(c_tmp*max(hmin,0.0)/max(hc,1.0e-06)+d_tmp, 1436 1296 min(c_tmp*max(hc,0.)/max(hmax,1.0e-06)+d_tmp, 1.0)) 1437 1297 ); 1438 //printf("%e, %e, \n", c_tmp, d_tmp);1439 //hfactor= max(0., min(5.0*max(hmin,0.0)/max(hmax,1.0e-06)-0.5, 1.0)1440 // );1441 //hfactor=1.0;1442 1298 // Set hfactor to zero smothly as hmin--> minimum_allowed_height. This 1443 1299 // avoids some 'chatter' for very shallow flows … … 1475 1331 1476 1332 beta_tmp = beta_w_dry + (beta_w - beta_w_dry) * hfactor; 1477 //beta_tmp = beta_w_dry*0. + (beta_w - beta_w_dry*0.) * hfactor;1478 //beta_tmp=1.0;1479 1480 1333 1481 1334 // Limit the gradient … … 1490 1343 // height 1491 1344 //----------------------------------- 1492 //hfactor=1.0;1493 //hfactor= max(0., min(5.0*max(hmin,0.0)/max(hc,1.0e-06)-0.5,1494 // min(5.0*max(hc,0.)/max(hmax,1.0e-06)-0.5, 1.0))1495 // );1496 1345 1497 1346 // Calculate the difference between vertex 0 of the auxiliary … … 1522 1371 1523 1372 // Limit the gradient 1373 // Same beta_tmp as for stage 1524 1374 //beta_tmp = beta_uh_dry + (beta_uh - beta_uh_dry) * hfactor; 1525 1375 limit_gradient(dqv, qmin, qmax, beta_tmp); … … 1531 1381 height_edge_values[k3+2] = height_centroid_values[k] + dqv[2]; 1532 1382 1533 1534 // REDEFINE hfactor for momentum terms -- make MORE first order1535 // Reduce beta linearly from 1-0 between depth ratio of 0.6-0.41536 //hfactor= max(0., min(5*max(hmin,0.0)/max(hc,1.0e-06)-2.0,1537 // min(5*max(hc,0.)/max(hmax,1.0e-06)-2.0, 1.0))1538 // );1539 //hfactor= max(0., min(5*max(hmin,0.0)/max(hmax,1.0e-06)-2.0,1540 // 1.0));1541 //hfactor=min( max(hmin,0.)/(max(hmin,0.)+10.*minimum_allowed_height), hfactor);1542 //hfactor=1.0;1543 1544 1383 //----------------------------------- 1545 1384 // xmomentum … … 1844 1683 if(extrapolate_velocity_second_order==1){ 1845 1684 //Convert velocity back to momenta at centroids 1846 xmom_centroid_values[k] = x mom_centroid_store[k];1847 ymom_centroid_values[k] = y mom_centroid_store[k];1685 xmom_centroid_values[k] = x_centroid_work[k]; 1686 ymom_centroid_values[k] = y_centroid_work[k]; 1848 1687 1849 1688 // Re-compute momenta at edges … … 1878 1717 1879 1718 } 1880 1881 free(xmom_centroid_store);1882 free(ymom_centroid_store);1883 //free(min_elevation_edgevalue);1884 //free(max_elevation_edgevalue);1885 //free(cell_wetness_scale);1886 //free(count_wet_neighbours);1887 1719 1888 1720 return 0; … … 1978 1810 *riverwall_elevation, 1979 1811 *riverwall_rowIndex, 1980 *riverwall_hydraulic_properties; 1812 *riverwall_hydraulic_properties, 1813 *edge_flux_work, 1814 *pressuregrad_work; 1981 1815 1982 1816 double timestep, epsilon, H0, g; … … 1984 1818 1985 1819 // Convert Python arguments to C 1986 if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOOOOiiOOOOOOOOiO ",1820 if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOOOOiiOOOOOOOOiOOO", 1987 1821 ×tep, 1988 1822 &epsilon, … … 2021 1855 &riverwall_rowIndex, 2022 1856 &ncol_riverwall_hydraulic_properties, 2023 &riverwall_hydraulic_properties 1857 &riverwall_hydraulic_properties, 1858 &edge_flux_work, 1859 &pressuregrad_work 2024 1860 )) { 2025 1861 report_python_error(AT, "could not parse input arguments"); … … 2059 1895 CHECK_C_CONTIG(riverwall_rowIndex); 2060 1896 CHECK_C_CONTIG(riverwall_hydraulic_properties); 1897 CHECK_C_CONTIG(edge_flux_work); 1898 CHECK_C_CONTIG(pressuregrad_work); 2061 1899 2062 1900 int number_of_elements = stage_edge_values -> dimensions[0]; … … 2103 1941 (long*) riverwall_rowIndex-> data, 2104 1942 ncol_riverwall_hydraulic_properties, 2105 (double*) riverwall_hydraulic_properties ->data); 1943 (double*) riverwall_hydraulic_properties ->data, 1944 (double*) edge_flux_work-> data, 1945 (double*) pressuregrad_work->data); 2106 1946 // Return updated flux timestep 2107 1947 return Py_BuildValue("d", timestep); … … 2274 2114 *ymom_vertex_values, 2275 2115 *elevation_vertex_values, 2276 *height_vertex_values; 2116 *height_vertex_values, 2117 *x_centroid_work, 2118 *y_centroid_work; 2277 2119 2278 2120 PyObject *domain; … … 2284 2126 2285 2127 // Convert Python arguments to C 2286 if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOOOOOOOOOii ",2128 if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOOOOOOOOOiiOO", 2287 2129 &domain, 2288 2130 &surrogate_neighbours, … … 2307 2149 &height_vertex_values, 2308 2150 &optimise_dry_cells, 2309 &extrapolate_velocity_second_order)) { 2151 &extrapolate_velocity_second_order, 2152 &x_centroid_work, 2153 &y_centroid_work)) { 2310 2154 2311 2155 report_python_error(AT, "could not parse input arguments"); … … 2334 2178 CHECK_C_CONTIG(elevation_vertex_values); 2335 2179 CHECK_C_CONTIG(height_vertex_values); 2180 CHECK_C_CONTIG(x_centroid_work); 2181 CHECK_C_CONTIG(y_centroid_work); 2336 2182 2337 2183 // Get the safety factor beta_w, set in the config.py file. … … 2384 2230 (double*) height_vertex_values -> data, 2385 2231 optimise_dry_cells, 2386 extrapolate_velocity_second_order); 2232 extrapolate_velocity_second_order, 2233 (double*) x_centroid_work -> data, 2234 (double*) y_centroid_work -> data); 2387 2235 2388 2236 if (e == -1) { -
trunk/anuga_work/development/gareth/tests/runup_sinusoid/runup_sinusoid.py
r9257 r9260 20 20 domain.set_name('runup_sinusoid_v2') # Output to file runup.sww 21 21 #domain.set_store_centroids(True) 22 domain.set_flow_algorithm('DE 0')22 domain.set_flow_algorithm('DE1') 23 23 domain.set_quantities_to_be_stored({'stage': 2, 'xmomentum': 2, 24 24 'ymomentum': 2, 'elevation': 2})
Note: See TracChangeset
for help on using the changeset viewer.
