Changeset 8992
- Timestamp:
- Sep 26, 2013, 3:24:18 PM (11 years ago)
- Location:
- trunk/anuga_work/development/gareth
- Files:
-
- 7 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/anuga_work/development/gareth/experimental/bal_and/swb2_domain.py
r8991 r8992 76 76 77 77 self.beta_w=1.0 78 self.beta_w_dry= 0.078 self.beta_w_dry=1.0 79 79 self.beta_uh=1.0 80 80 self.beta_uh_dry=0.0 -
trunk/anuga_work/development/gareth/experimental/bal_and/swb2_domain_ext.c
r8991 r8992 355 355 double speed_max_last, vol, smooth; 356 356 357 double *count_wet_edges, *edgeflux_store, *pressuregrad_store; 358 359 //count_wet_neighbours = malloc(number_of_elements*sizeof(double)); 360 count_wet_edges = malloc(number_of_elements*sizeof(double)); 357 double *edgeflux_store, *pressuregrad_store; 358 361 359 edgeflux_store = malloc(number_of_elements*9*sizeof(double)); 362 360 pressuregrad_store = malloc(number_of_elements*3*sizeof(double)); … … 445 443 446 444 447 //if((count_wet_edges[k]<=1.0 | count_wet_neighbours[k]<=1.0) &448 // (n>=0 && (count_wet_edges[n]<=1.0 | count_wet_neighbours[n]<=1.0))){449 // smooth=0.0;450 //}451 445 452 446 // Edge flux computation (triangle k, edge i) … … 474 468 pressuregrad_store[ki]=bedslope_work; 475 469 476 //if((k== 2 | n==2) | (k==4|n==4)){477 // if(n<0){470 //if((k==5 | n==5) ){ 471 ////if(n<0){ 478 472 // printf("k n: %d,%d \n", k, n); 479 473 // printf("z's: %e,%e \n", zl, zr); 480 474 // printf("h_LR's: %e,%e \n", h_left,h_right); 481 475 // printf("h's: %e,%e \n", hle,hre); 476 // printf("hc's: %e,%e \n", hc,hc_n); 477 // printf("qc's: %e,%e, %e, %e \n", ql[1], qr[1],ql[2], qr[2]); 482 478 // printf("pf: %e \n", pressure_flux); 483 479 // printf("flux: %e, %e, %e\n", edgeflux[0], edgeflux[1], edgeflux[2]); … … 602 598 // GD HACK 603 599 // Option to limit advective fluxes 604 if(hc > -h0){600 //if(hc > -h0){ 605 601 stage_explicit_update[k] += edgeflux_store[ki3+0]; 606 602 // Cut these terms for shallow flows, and protect stationary states from round-off error 607 603 xmom_explicit_update[k] += edgeflux_store[ki3+1]; 608 604 ymom_explicit_update[k] += edgeflux_store[ki3+2]; 609 }else{610 stage_explicit_update[k] += edgeflux_store[ki3+0];611 }605 //}else{ 606 // stage_explicit_update[k] += edgeflux_store[ki3+0]; 607 //} 612 608 613 609 … … 619 615 // GD HACK 620 616 // Compute bed slope term 621 if(hc > -h0){617 //if(hc > -h0){ 622 618 xmom_explicit_update[k] -= normals[ki2]*pressuregrad_store[ki]; 623 619 ymom_explicit_update[k] -= normals[ki2+1]*pressuregrad_store[ki]; 624 }else{625 xmom_explicit_update[k] *= 0.;626 ymom_explicit_update[k] *= 0.;627 }620 //}else{ 621 // xmom_explicit_update[k] *= 0.; 622 // ymom_explicit_update[k] *= 0.; 623 //} 628 624 629 625 } // end edge i … … 694 690 //return -1; 695 691 // 696 free(count_wet_edges);697 692 //free(count_wet_neighbours); 698 693 free(edgeflux_store); … … 904 899 //count_wet_neighbours = malloc(number_of_elements*sizeof(int)); 905 900 906 if(extrapolate_velocity_second_order==1){907 // Replace momentum centroid with velocity centroid to allow velocity908 // extrapolation This will be changed back at the end of the routine909 for (k=0; k<number_of_elements; k++){910 911 dk = max(stage_centroid_values[k]-elevation_centroid_values[k],minimum_allowed_height);912 xmom_centroid_store[k] = xmom_centroid_values[k];913 xmom_centroid_values[k] = xmom_centroid_values[k]/dk;914 915 ymom_centroid_store[k] = ymom_centroid_values[k];916 ymom_centroid_values[k] = ymom_centroid_values[k]/dk;917 918 min_elevation_edgevalue[k] = min(elevation_edge_values[3*k],919 min(elevation_edge_values[3*k+1],920 elevation_edge_values[3*k+2]));921 max_elevation_edgevalue[k] = max(elevation_edge_values[3*k],922 max(elevation_edge_values[3*k+1],923 elevation_edge_values[3*k+2]));924 925 // SET HEIGHT IN PLACE926 height_centroid_values[k] = max(stage_centroid_values[k] - elevation_centroid_values[k], 0.);927 928 }929 930 }931 932 901 // 933 902 // Compute some useful statistics on wetness/dryness … … 937 906 // Check if cell k is wet 938 907 //if(stage_centroid_values[k] > elevation_centroid_values[k]){ 939 if(stage_centroid_values[k] > elevation_centroid_values[k] + 2 *minimum_allowed_height+epsilon){908 if(stage_centroid_values[k] > elevation_centroid_values[k] + 2.*minimum_allowed_height+epsilon){ 940 909 //if(stage_centroid_values[k] > max_elevation_edgevalue[k] + minimum_allowed_height+epsilon){ 941 910 cell_wetness_scale[k] = 1.; 942 911 } 943 912 } 913 914 // Alternative 'PROTECT' step 915 for(k=0; k<number_of_elements;k++){ 916 if((cell_wetness_scale[k]==0. ) ){ 917 xmom_centroid_values[k]=0.; 918 ymom_centroid_values[k]=0.; 919 } 920 } 921 944 922 // 945 923 //for(k=0; k<number_of_elements;k++){ … … 959 937 //} 960 938 961 // Alternative 'PROTECT' step 962 for(k=0; k<number_of_elements;k++){ 963 if(1==1){ 964 if((cell_wetness_scale[k]==0. ) ){ 965 xmom_centroid_store[k]=0.; 966 xmom_centroid_values[k]=0.; 967 ymom_centroid_store[k]=0.; 968 ymom_centroid_values[k]=0.; 969 } 970 } 971 } 972 939 if(extrapolate_velocity_second_order==1){ 940 // Replace momentum centroid with velocity centroid to allow velocity 941 // extrapolation This will be changed back at the end of the routine 942 for (k=0; k<number_of_elements; k++){ 943 944 dk = max(stage_centroid_values[k]-elevation_centroid_values[k],minimum_allowed_height); 945 xmom_centroid_store[k] = xmom_centroid_values[k]; 946 xmom_centroid_values[k] = xmom_centroid_values[k]/dk; 947 948 ymom_centroid_store[k] = ymom_centroid_values[k]; 949 ymom_centroid_values[k] = ymom_centroid_values[k]/dk; 950 951 min_elevation_edgevalue[k] = min(elevation_edge_values[3*k], 952 min(elevation_edge_values[3*k+1], 953 elevation_edge_values[3*k+2])); 954 max_elevation_edgevalue[k] = max(elevation_edge_values[3*k], 955 max(elevation_edge_values[3*k+1], 956 elevation_edge_values[3*k+2])); 957 958 // SET HEIGHT IN PLACE 959 height_centroid_values[k] = max(stage_centroid_values[k] - elevation_centroid_values[k], 0.); 960 961 } 962 963 } 964 973 965 // First treat all 'fully wet' cells, then treat 'partially dry' cells 974 966 //for(k_wetdry=0; k_wetdry<2; k_wetdry++){ … … 1076 1068 //if(stage_centroid_values[k2]<=max_elevation_edgevalue[k2]+epsilon){ 1077 1069 //if(stage_centroid_values[k2]<=elevation_centroid_values[k2]+minimum_allowed_height+epsilon){ 1078 if(k2<0 || cell_wetness_scale[k2]==0.0){ 1070 // FIXME: Remove cell_wetness_scale if you don't need it 1071 if(k2<0 || (cell_wetness_scale[k2]==-10.0 && stage_centroid_values[k]<stage_centroid_values[k2]) ){ 1079 1072 k2 = k ; 1080 1073 } 1081 1074 //if(stage_centroid_values[k0]<=max_elevation_edgevalue[k0]+epsilon){ 1082 1075 //if(stage_centroid_values[k0]<=elevation_centroid_values[k0]+minimum_allowed_height+epsilon){ 1083 if(k0 < 0 || cell_wetness_scale[k0]==0.0){1076 if(k0 < 0 || (cell_wetness_scale[k0]==-10.0 && stage_centroid_values[k]<stage_centroid_values[k0])){ 1084 1077 k0 = k ; 1085 1078 } 1086 1079 //if(stage_centroid_values[k1]<=max_elevation_edgevalue[k1]+epsilon){ 1087 1080 //if(stage_centroid_values[k1]<=elevation_centroid_values[k1]+minimum_allowed_height+epsilon){ 1088 if(k1 < 0 || cell_wetness_scale[k1]==0.0){1081 if(k1 < 0 || (cell_wetness_scale[k1]==-10.0 && stage_centroid_values[k]<stage_centroid_values[k1])){ 1089 1082 k1 = k ; 1090 1083 } … … 1135 1128 1136 1129 // Isolated wet cell -- constant stage/depth extrapolation 1137 stage_edge_values[k3] = stage_centroid_values[k];1138 stage_edge_values[k3+1] = stage_centroid_values[k];1139 stage_edge_values[k3+2] = stage_centroid_values[k];1130 //stage_edge_values[k3] = stage_centroid_values[k]; 1131 //stage_edge_values[k3+1] = stage_centroid_values[k]; 1132 //stage_edge_values[k3+2] = stage_centroid_values[k]; 1140 1133 1141 1134 dk=max(stage_centroid_values[k]-elevation_centroid_values[k],0.); … … 1143 1136 height_edge_values[k3+1] = dk; 1144 1137 height_edge_values[k3+2] = dk; 1138 1139 stage_edge_values[k3] = elevation_centroid_values[k]+dk; 1140 stage_edge_values[k3+1] = elevation_centroid_values[k]+dk; 1141 stage_edge_values[k3+2] = elevation_centroid_values[k]+dk; 1145 1142 1146 1143 xmom_edge_values[k3] = xmom_centroid_values[k]; … … 1614 1611 // Re-compute momenta at edges 1615 1612 for (i=0; i<3; i++){ 1616 de[i] = height_edge_values[k3+i]; //max(stage_edge_values[k3+i]-elevation_edge_values[k3+i],0.0); 1613 //de[i] = height_edge_values[k3+i]; //max(stage_vertex_values[k3+i]-elevation_vertex_values[k3+i],0.0); 1614 // Method to perhaps slightly reduce momenta 1615 de[i] = min(height_edge_values[k3+i], max(stage_edge_values[k3+i]-elevation_edge_values[k3+i],0.)); //max(stage_edge_values[k3+i]-elevation_edge_values[k3+i],0.0); 1617 1616 //if(de[i]<minimum_allowed_height){ 1618 1617 // de[i]=0.; … … 1624 1623 // Re-compute momenta at vertices 1625 1624 for (i=0; i<3; i++){ 1626 de[i] = height_vertex_values[k3+i]; //max(stage_vertex_values[k3+i]-elevation_vertex_values[k3+i],0.0); 1625 //de[i] = height_vertex_values[k3+i]; //max(stage_vertex_values[k3+i]-elevation_vertex_values[k3+i],0.0); 1626 // Method to perhaps slightly reduce momenta 1627 de[i] = min(height_vertex_values[k3+i], max(stage_vertex_values[k3+i]-elevation_vertex_values[k3+i],0.)); //max(stage_vertex_values[k3+i]-elevation_vertex_values[k3+i],0.0); 1627 1628 xmom_vertex_values[k3+i]=xmom_vertex_values[k3+i]*de[i]; 1628 1629 ymom_vertex_values[k3+i]=ymom_vertex_values[k3+i]*de[i]; -
trunk/anuga_work/development/gareth/tests/parabolic/parabolic.py
r8991 r8992 23 23 domain.set_name('parabola_v2') # Output to file runup.sww 24 24 domain.set_datadir('.') # Use current folder 25 domain.set_minimum_allowed_height(0.01)25 #domain.set_minimum_allowed_height(0.01) 26 26 27 27 #domain.set_flow_algorithm('tsunami') -
trunk/anuga_work/development/gareth/tests/runup/runuplot.py
r8353 r8992 8 8 import scipy 9 9 from matplotlib import pyplot as pyplot 10 import util # Routines to read in and work with ANUGA output 10 #import util # Routines to read in and work with ANUGA output 11 from bal_and import plot_utils as util 11 12 12 13 p2=util.get_output('runup_v2.sww', minimum_allowed_height=1.0e-03) … … 58 59 pyplot.plot(p.x[v],p.stage[5,v]) 59 60 pyplot.plot(p.x[v],p.stage[5,v],'o') 60 pyplot.plot(p.x[v],p.elev[ v])61 pyplot.plot(p.x[v],p.elev[5,v]) 61 62 pyplot.xlabel('x (m)') 62 63 pyplot.ylabel('z (m)') … … 91 92 92 93 pyplot.clf() 93 pyplot.scatter(p.x,p.y,c=p.elev ,edgecolors='none', s=25)94 pyplot.scatter(p.x,p.y,c=p.elev[1,],edgecolors='none', s=25) 94 95 pyplot.colorbar() 95 96 #pyplot.quiver(x_cent,y_cent,xvel_cent[15,:],yvel_cent[15,:]) … … 100 101 101 102 pyplot.clf() 102 pyplot.scatter(p.x,p.y,c=p.elev ,edgecolors='none', s=25)103 pyplot.scatter(p.x,p.y,c=p.elev[1,],edgecolors='none', s=25) 103 104 pyplot.colorbar() 104 105 #pyplot.quiver(x_cent,y_cent,xvel_cent[150,:],yvel_cent[150,:]) … … 111 112 pyplot.plot(p.x[v],p.stage[150,v]) 112 113 pyplot.plot(p.x[v],p.stage[150,v],'o') 113 pyplot.plot(p.x[v],p.elev[ v])114 pyplot.plot(p.x[v],p.elev[150,v]) 114 115 pyplot.xlabel('x (m)') 115 116 pyplot.ylabel('z (m)') -
trunk/anuga_work/development/gareth/tests/runup_sinusoid/runup_sinusoid.py
r8991 r8992 44 44 init_ws=-0.2 45 45 bumpiness=50. # Higher = shorter wavelength oscillations in topography 46 tstep=0. 0247 lasttime= 40.46 tstep=0.2 47 lasttime=90. 48 48 49 49 #domain.minimum_allowed_height=domain.minimum_allowed_height*scale_me # Seems needed to make the algorithms behave -
trunk/anuga_work/development/gareth/tests/runup_sinusoid/runup_sinusoidplot.py
r8990 r8992 45 45 t2=len(p1.time)-1 46 46 47 pyplot.scatter(p1.x,p1.y,c=p1.elev ,edgecolors='none', s=25)47 pyplot.scatter(p1.x,p1.y,c=p1.elev[1,:],edgecolors='none', s=25) 48 48 pyplot.colorbar() 49 49 pyplot.quiver(p1.x,p1.y,p1.xvel[t1,:],p1.yvel[t1,:]) … … 53 53 pyplot.clf() 54 54 # Plot vertex values 55 pyplot.scatter(p2.x,p2.y,c=p2.elev ,edgecolors='none', s=25)55 pyplot.scatter(p2.x,p2.y,c=p2.elev[1,:],edgecolors='none', s=25) 56 56 pyplot.colorbar() 57 57 pyplot.quiver(p2.x,p2.y,p2.xvel[t1,:],p2.yvel[t1,:]) … … 61 61 pyplot.clf() 62 62 # Plot vertex values 63 pyplot.scatter(p1.x,p1.y,c=p1.elev ,edgecolors='none', s=25)63 pyplot.scatter(p1.x,p1.y,c=p1.elev[1,:],edgecolors='none', s=25) 64 64 pyplot.colorbar() 65 65 pyplot.quiver(p1.x,p1.y,p1.xvel[t2,:],p1.yvel[t2,:]) … … 69 69 pyplot.clf() 70 70 # Plot vertex values 71 pyplot.scatter(p2.x,p2.y,c=p2.elev ,edgecolors='none', s=25)71 pyplot.scatter(p2.x,p2.y,c=p2.elev[1,:],edgecolors='none', s=25) 72 72 pyplot.colorbar() 73 73 pyplot.quiver(p2.x,p2.y,p2.xvel[t2,:],p2.yvel[t2,:]) -
trunk/anuga_work/development/gareth/tests/shallow_steep_slope/channel_SU_sparse.py
r8990 r8992 43 43 return stg#*(stg>topo) + topo*(stg<=topo) 44 44 45 line1=[ [ 9.,0.], [9., 100.] ]45 line1=[ [0.,0.], [0., 100.] ] 46 46 Qin=0.5 47 47 Inlet_operator(domain, line1,Qin)
Note: See TracChangeset
for help on using the changeset viewer.