Changeset 8992


Ignore:
Timestamp:
Sep 26, 2013, 3:24:18 PM (11 years ago)
Author:
davies
Message:

Audusse improvements

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  
    7676
    7777        self.beta_w=1.0
    78         self.beta_w_dry=0.0
     78        self.beta_w_dry=1.0
    7979        self.beta_uh=1.0
    8080        self.beta_uh_dry=0.0
  • trunk/anuga_work/development/gareth/experimental/bal_and/swb2_domain_ext.c

    r8991 r8992  
    355355    double speed_max_last, vol, smooth;
    356356
    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
    361359    edgeflux_store = malloc(number_of_elements*9*sizeof(double));
    362360    pressuregrad_store = malloc(number_of_elements*3*sizeof(double));
     
    445443           
    446444
    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             //}
    451445           
    452446            // Edge flux computation (triangle k, edge i)
     
    474468            pressuregrad_store[ki]=bedslope_work;
    475469           
    476             //if((k==2 | n==2) | (k==4|n==4)){
    477             //if(n<0){
     470            //if((k==5 | n==5) ){
     471            ////if(n<0){
    478472            //    printf("k n: %d,%d \n", k, n);
    479473            //    printf("z's: %e,%e \n", zl, zr);
    480474            //    printf("h_LR's: %e,%e \n", h_left,h_right);
    481475            //    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]);
    482478            //    printf("pf: %e \n", pressure_flux);
    483479            //    printf("flux: %e, %e, %e\n", edgeflux[0], edgeflux[1], edgeflux[2]);
     
    602598            // GD HACK
    603599            // Option to limit advective fluxes
    604             if(hc > -h0){
     600            //if(hc > -h0){
    605601                stage_explicit_update[k] += edgeflux_store[ki3+0];
    606602                // Cut these terms for shallow flows, and protect stationary states from round-off error
    607603                xmom_explicit_update[k] += edgeflux_store[ki3+1];
    608604                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            //}
    612608
    613609
     
    619615            // GD HACK
    620616            // Compute bed slope term
    621             if(hc > -h0){
     617            //if(hc > -h0){
    622618                xmom_explicit_update[k] -= normals[ki2]*pressuregrad_store[ki];
    623619                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            //}
    628624
    629625        } // end edge i
     
    694690    //return -1;
    695691//
    696     free(count_wet_edges);
    697692    //free(count_wet_neighbours);
    698693    free(edgeflux_store);
     
    904899  //count_wet_neighbours = malloc(number_of_elements*sizeof(int));
    905900 
    906   if(extrapolate_velocity_second_order==1){
    907       // Replace momentum centroid with velocity centroid to allow velocity
    908       // extrapolation This will be changed back at the end of the routine
    909       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 PLACE
    926           height_centroid_values[k] = max(stage_centroid_values[k] - elevation_centroid_values[k], 0.);
    927            
    928           }
    929 
    930       }
    931 
    932901  //
    933902  // Compute some useful statistics on wetness/dryness
     
    937906      // Check if cell k is wet
    938907      //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){
    940909      //if(stage_centroid_values[k] > max_elevation_edgevalue[k] + minimum_allowed_height+epsilon){
    941910          cell_wetness_scale[k] = 1.; 
    942911      }
    943912  }
     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 
    944922  //
    945923  //for(k=0; k<number_of_elements;k++){
     
    959937  //}
    960938
    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
    973965  // First treat all 'fully wet' cells, then treat 'partially dry' cells
    974966  //for(k_wetdry=0; k_wetdry<2; k_wetdry++){
     
    10761068          //if(stage_centroid_values[k2]<=max_elevation_edgevalue[k2]+epsilon){
    10771069          //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]) ){
    10791072              k2 = k ;
    10801073          }
    10811074          //if(stage_centroid_values[k0]<=max_elevation_edgevalue[k0]+epsilon){
    10821075          //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])){
    10841077              k0 = k ;
    10851078          }
    10861079          //if(stage_centroid_values[k1]<=max_elevation_edgevalue[k1]+epsilon){
    10871080          //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])){
    10891082              k1 = k ;
    10901083          }
     
    11351128
    11361129              // 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];
    11401133
    11411134              dk=max(stage_centroid_values[k]-elevation_centroid_values[k],0.);
     
    11431136              height_edge_values[k3+1] = dk;
    11441137              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;
    11451142
    11461143              xmom_edge_values[k3]    = xmom_centroid_values[k];
     
    16141611          // Re-compute momenta at edges
    16151612          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);
    16171616              //if(de[i]<minimum_allowed_height){
    16181617              //  de[i]=0.;
     
    16241623          // Re-compute momenta at vertices
    16251624          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);
    16271628              xmom_vertex_values[k3+i]=xmom_vertex_values[k3+i]*de[i];
    16281629              ymom_vertex_values[k3+i]=ymom_vertex_values[k3+i]*de[i];
  • trunk/anuga_work/development/gareth/tests/parabolic/parabolic.py

    r8991 r8992  
    2323domain.set_name('parabola_v2')                         # Output to file runup.sww
    2424domain.set_datadir('.')                          # Use current folder
    25 domain.set_minimum_allowed_height(0.01)
     25#domain.set_minimum_allowed_height(0.01)
    2626
    2727#domain.set_flow_algorithm('tsunami')
  • trunk/anuga_work/development/gareth/tests/runup/runuplot.py

    r8353 r8992  
    88import scipy
    99from 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
     11from bal_and import plot_utils as util
    1112
    1213p2=util.get_output('runup_v2.sww', minimum_allowed_height=1.0e-03)
     
    5859pyplot.plot(p.x[v],p.stage[5,v])
    5960pyplot.plot(p.x[v],p.stage[5,v],'o')
    60 pyplot.plot(p.x[v],p.elev[v])
     61pyplot.plot(p.x[v],p.elev[5,v])
    6162pyplot.xlabel('x (m)')
    6263pyplot.ylabel('z (m)')
     
    9192
    9293pyplot.clf()
    93 pyplot.scatter(p.x,p.y,c=p.elev,edgecolors='none', s=25)
     94pyplot.scatter(p.x,p.y,c=p.elev[1,],edgecolors='none', s=25)
    9495pyplot.colorbar()
    9596#pyplot.quiver(x_cent,y_cent,xvel_cent[15,:],yvel_cent[15,:])
     
    100101
    101102pyplot.clf()
    102 pyplot.scatter(p.x,p.y,c=p.elev,edgecolors='none', s=25)
     103pyplot.scatter(p.x,p.y,c=p.elev[1,],edgecolors='none', s=25)
    103104pyplot.colorbar()
    104105#pyplot.quiver(x_cent,y_cent,xvel_cent[150,:],yvel_cent[150,:])
     
    111112pyplot.plot(p.x[v],p.stage[150,v])
    112113pyplot.plot(p.x[v],p.stage[150,v],'o')
    113 pyplot.plot(p.x[v],p.elev[v])
     114pyplot.plot(p.x[v],p.elev[150,v])
    114115pyplot.xlabel('x (m)')
    115116pyplot.ylabel('z (m)')
  • trunk/anuga_work/development/gareth/tests/runup_sinusoid/runup_sinusoid.py

    r8991 r8992  
    4444init_ws=-0.2
    4545bumpiness=50. # Higher = shorter wavelength oscillations in topography
    46 tstep=0.02
    47 lasttime=40.
     46tstep=0.2
     47lasttime=90.
    4848
    4949#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  
    4545t2=len(p1.time)-1
    4646
    47 pyplot.scatter(p1.x,p1.y,c=p1.elev,edgecolors='none', s=25)
     47pyplot.scatter(p1.x,p1.y,c=p1.elev[1,:],edgecolors='none', s=25)
    4848pyplot.colorbar()
    4949pyplot.quiver(p1.x,p1.y,p1.xvel[t1,:],p1.yvel[t1,:])
     
    5353pyplot.clf()
    5454# Plot vertex values
    55 pyplot.scatter(p2.x,p2.y,c=p2.elev,edgecolors='none', s=25)
     55pyplot.scatter(p2.x,p2.y,c=p2.elev[1,:],edgecolors='none', s=25)
    5656pyplot.colorbar()
    5757pyplot.quiver(p2.x,p2.y,p2.xvel[t1,:],p2.yvel[t1,:])
     
    6161pyplot.clf()
    6262# Plot vertex values
    63 pyplot.scatter(p1.x,p1.y,c=p1.elev,edgecolors='none', s=25)
     63pyplot.scatter(p1.x,p1.y,c=p1.elev[1,:],edgecolors='none', s=25)
    6464pyplot.colorbar()
    6565pyplot.quiver(p1.x,p1.y,p1.xvel[t2,:],p1.yvel[t2,:])
     
    6969pyplot.clf()
    7070# Plot vertex values
    71 pyplot.scatter(p2.x,p2.y,c=p2.elev,edgecolors='none', s=25)
     71pyplot.scatter(p2.x,p2.y,c=p2.elev[1,:],edgecolors='none', s=25)
    7272pyplot.colorbar()
    7373pyplot.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  
    4343    return stg#*(stg>topo) + topo*(stg<=topo)
    4444
    45 line1=[ [9.,0.], [9., 100.] ]
     45line1=[ [0.,0.], [0., 100.] ]
    4646Qin=0.5
    4747Inlet_operator(domain, line1,Qin)
Note: See TracChangeset for help on using the changeset viewer.