Changeset 8549


Ignore:
Timestamp:
Sep 1, 2012, 7:44:08 PM (13 years ago)
Author:
davies
Message:

Bugfixes for balanced_dev

Location:
trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/anuga_core/source/anuga/structures/inlet_operator.py

    r8506 r8549  
    6262        volume = Q*timestep
    6363       
     64        #print Q, volume
     65       
    6466        assert current_volume + volume >= 0.0, 'Requesting too much water to be removed from an inlet!'
    6567

  • trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain_ext.c

    r8547 r8549  
    196196  // Flux formulas
    197197  flux_left[0] = u_left*h_left;
    198   flux_left[1] = u_left*uh_left + 0.0*g*h_left*h_left;
     198  flux_left[1] = u_left*uh_left ; //+ 0.5*g*h_left*h_left;
    199199  flux_left[2] = u_left*vh_left;
    200200
    201201  flux_right[0] = u_right*h_right;
    202   flux_right[1] = u_right*uh_right + 0.0*g*h_right*h_right;
     202  flux_right[1] = u_right*uh_right ; //+ 0.5*g*h_right*h_right;
    203203  flux_right[2] = u_right*vh_right;
    204204
    … …  
    221221      edgeflux[i] *= inverse_denominator;
    222222    }
    223 
     223    // Separate pressure flux, so we can apply different wet-dry hacks to it
    224224    *pressure_flux = 0.5*g*( s_max*h_left*h_left -s_min*h_right*h_right)*inverse_denominator;
    225225   
    … …  
    288288
    289289    max_bed_edgevalue = malloc(number_of_elements*sizeof(double));
    290     min_bed_edgevalue = malloc(number_of_elements*sizeof(double));
     290    //min_bed_edgevalue = malloc(number_of_elements*sizeof(double));
    291291    // Start computation
    292292    call++; // Flag 'id' of flux calculation for this timestep
    … …  
    304304                                   max(bed_edge_values[3*k+1], bed_edge_values[3*k+2]));
    305305        //max_bed_edgevalue[k]= 0.5*(max_bed_edgevalue[j]+bed_centroid_values[k]);
    306         min_bed_edgevalue[k] = min(bed_edge_values[3*k],
    307                                    min(bed_edge_values[3*k+1], bed_edge_values[3*k+2]));
     306        //min_bed_edgevalue[k] = min(bed_edge_values[3*k],
     307        //                           min(bed_edge_values[3*k+1], bed_edge_values[3*k+2]));
    308308   
    309309    }
    … …  
    408408            // Prevent outflow from 'seriously' dry cells
    409409            // Idea: The cell will not go dry if:
    410             // mass_flux = edgeflux[0]*(dt*edgelength) <= Area_triangle*hc
     410            // mass_flux = edgeflux[0]*(dt*edgelength) <= (1/3)Area_triangle*hc
    411411            if((stage_centroid_values[k]<=max_bed_edgevalue[k])|
    412412               (ql[0]<=zl)){
    … …  
    469469                ymom_explicit_update[k] -= normals[ki2+1]*bedslope_work;
    470470            }else{
    471                 bedslope_work = 0.;//-pressure_flux*length ;//g*length*( 0.*(ql[0])-0.5*max(ql[0]-zl,0.)*(ql[0]-zl) );
     471                //bedslope_work = 0.;//-pressure_flux*length ;//g*length*( 0.*(ql[0])-0.5*max(ql[0]-zl,0.)*(ql[0]-zl) );
    472472                xmom_explicit_update[k] *= 0.;
    473473                ymom_explicit_update[k] *= 0.;
    … …  
    524524                    ymom_explicit_update[n] += normals[ki2+1]*bedslope_work;
    525525                }else{
    526                     bedslope_work = 0.;//-pressure_flux*length;//g*length*(0.*(qr[0])-0.5*max(qr[0]-zr,0.)*(qr[0]-zr));
     526                    //bedslope_work = 0.;//-pressure_flux*length;//g*length*(0.*(qr[0])-0.5*max(qr[0]-zr,0.)*(qr[0]-zr));
    527527                    xmom_explicit_update[n] *= 0.;
    528528                    ymom_explicit_update[n] *= 0.;
    … …  
    633633//
    634634    free(max_bed_edgevalue);
    635     free(min_bed_edgevalue);
     635    //free(min_bed_edgevalue);
    636636
    637637    return timestep;
    … …  
    671671
    672672      if (hc < max(minimum_relative_height*(bmax-zc[k]), minimum_allowed_height)) {
    673        
    674673        // Set momentum to zero and ensure h is non negative
    675         // NOTE: THIS IS IMPORTANT -- WE ARE SETTING MOMENTUM TO ZERO
    676         //if(hc<=epsilon){
    677674            xmomc[k] = 0.0;
    678675            ymomc[k] = 0.0;
    679         //}
    680676
    681677        if (hc <= 0.0){
    682678             // Definine the minimum bed edge value on triangle k.
    683              // Setting = minimum edge value can lead to mass conservation problems
    684679             //bmin = 0.5*min((zv[3*k]+zv[3*k+1]),min((zv[3*k+1]+zv[3*k+2]),(zv[3*k+2]+zv[3*k])));
    685680             //bmin =0.5*bmin + 0.5*min(zv[3*k],min(zv[3*k+1],zv[3*k+2]));
    686              // Setting = minimum vertex value seems better, but might tend to be less smooth
    687681             //bmin =min(zv[3*k],min(zv[3*k+1],zv[3*k+2])) -minimum_allowed_height;
    688682             bmin=zc[k]-minimum_allowed_height;
    689683             // Minimum allowed stage = bmin
     684
    690685             // WARNING: ADDING MASS if wc[k]<bmin
    691686             if(wc[k]<bmin){
    … …  
    697692             
    698693
    699                  // Set vertex values as well. Seems that this shouldn't be
    700                  // needed. However from memory this is important at the first
     694                 // FIXME: Set vertex values as well. Seems that this shouldn't be
     695                 // needed. However, from memory this is important at the first
    701696                 // time step, for 'dry' areas where the designated stage is
    702697                 // less than the bed centroid value
    … …  
    813808 
    814809  double *xmom_centroid_store, *ymom_centroid_store, *stage_centroid_store, *min_elevation_edgevalue, *max_elevation_edgevalue;
    815   double *neighbour_wetness_scale;
     810  double *cell_wetness_scale;
    816811  int *count_wet_neighbours;
    817812
    … …  
    821816  ymom_centroid_store = malloc(number_of_elements*sizeof(double));
    822817  stage_centroid_store = malloc(number_of_elements*sizeof(double));
    823   min_elevation_edgevalue = malloc(number_of_elements*sizeof(double));
     818  //min_elevation_edgevalue = malloc(number_of_elements*sizeof(double));
    824819  max_elevation_edgevalue = malloc(number_of_elements*sizeof(double));
    825   neighbour_wetness_scale = malloc(number_of_elements*sizeof(double));
     820  cell_wetness_scale = malloc(number_of_elements*sizeof(double));
    826821  count_wet_neighbours = malloc(number_of_elements*sizeof(int));
    827822 
    … …  
    838833          ymom_centroid_values[k] = ymom_centroid_values[k]/dk;
    839834
    840           min_elevation_edgevalue[k] = min(elevation_edge_values[3*k],
    841                                            min(elevation_edge_values[3*k+1],
    842                                                elevation_edge_values[3*k+2]));
     835          //min_elevation_edgevalue[k] = min(elevation_edge_values[3*k],
     836          //                                 min(elevation_edge_values[3*k+1],
     837          //                                     elevation_edge_values[3*k+2]));
    843838          max_elevation_edgevalue[k] = max(elevation_edge_values[3*k],
    844839                                           max(elevation_edge_values[3*k+1],
    … …  
    852847  for(k=0; k<number_of_elements;k++){
    853848      count_wet_neighbours[k]=0;
    854       neighbour_wetness_scale[k] = 0.;
    855       // Cell k is wet
    856       if(stage_centroid_values[k] > max_elevation_edgevalue[k] + minimum_allowed_height+epsilon){
    857           neighbour_wetness_scale[k] = 1.; 
     849      cell_wetness_scale[k] = 0.;
     850      // Check if cell k is wet
     851      if(stage_centroid_values[k] > elevation_centroid_values[k] + minimum_allowed_height+epsilon){
     852          cell_wetness_scale[k] = 1.; 
    858853      }
    859854      // Count the wet neighbours
    … …  
    871866  for(k_wetdry=0; k_wetdry<2; k_wetdry++){
    872867      //if(((stage_centroid_values[k] > max_elevation_edgevalue[k]))==k_wetdry){
    873       if(((stage_centroid_values[k] < elevation_centroid_values[k]+minimum_allowed_height))==k_wetdry){
    874         continue;
    875       }
     868      //if(((stage_centroid_values[k] < elevation_centroid_values[k]+minimum_allowed_height))==k_wetdry){
    876869      // For partially wet cells, we now know that the edges of neighbouring
    877870      // fully wet cells have been defined
    … …  
    880873      for (k = 0; k < number_of_elements; k++)
    881874      {
     875
     876        // First treat all 'fully wet' cells, then treat 'partially dry' cells
     877        // For partially wet cells, we now know that the edges of neighbouring
     878        // fully wet cells have been defined
     879        if( cell_wetness_scale[k]==(1.0*k_wetdry) ){
     880          continue;
     881        }
     882
     883        // Useful indices
    882884        k3=k*3;
    883885        k6=k*6;
    … …  
    950952          // used
    951953          //if(stage_centroid_values[k2]<=max_elevation_edgevalue[k2]+epsilon){
    952           if(stage_centroid_values[k2]<=elevation_centroid_values[k2]+minimum_allowed_height+epsilon){
     954          //if(stage_centroid_values[k2]<=elevation_centroid_values[k2]+minimum_allowed_height+epsilon){
     955          if(cell_wetness_scale[k2]==0.0){
    953956              k2 = k ;
    954957          }
    955958          //if(stage_centroid_values[k0]<=max_elevation_edgevalue[k0]+epsilon){
    956           if(stage_centroid_values[k0]<=elevation_centroid_values[k0]+minimum_allowed_height+epsilon){
     959          //if(stage_centroid_values[k0]<=elevation_centroid_values[k0]+minimum_allowed_height+epsilon){
     960          if(cell_wetness_scale[k0]==0.0){
    957961              k0 = k ;
    958962          }
    959963          //if(stage_centroid_values[k1]<=max_elevation_edgevalue[k1]+epsilon){
    960           if(stage_centroid_values[k1]<=elevation_centroid_values[k1]+minimum_allowed_height+epsilon){
     964          //if(stage_centroid_values[k1]<=elevation_centroid_values[k1]+minimum_allowed_height+epsilon){
     965          if(cell_wetness_scale[k1]==0.0){
    961966              k1 = k ;
    962967          }
    … …  
    9991004           
    10001005          // Treat triangles with zero or 1 wet neighbours.
    1001           if ((area2 <= 0)) //|(count_wet_neighbours[k]==0))
     1006          if ((area2 <= 0)|(cell_wetness_scale[k]==0.)) //|(count_wet_neighbours[k]==0))
    10021007          {
    1003               if( ((k==k0) & (k==k1) & (k==k2))){
    1004                   // Isolated wet cell  -- use constant depth extrapolation for stage
    1005                   //stage_edge_values[k3]   = stage_centroid_values[k]- elevation_centroid_values[k] + elevation_edge_values[k3];
    1006                   //stage_edge_values[k3+1] = stage_centroid_values[k]- elevation_centroid_values[k] + elevation_edge_values[k3+1];
    1007                   //stage_edge_values[k3+2] = stage_centroid_values[k]- elevation_centroid_values[k] + elevation_edge_values[k3+2];
    1008 
    1009                   // Isolated wet cell -- constant depth extrapolation, but don't let water flow uphill
    1010                   //stage_edge_values[k3]   = min(stage_centroid_values[k], stage_centroid_values[k]- elevation_centroid_values[k] + elevation_edge_values[k3]);
    1011                   //stage_edge_values[k3+1] = min(stage_centroid_values[k], stage_centroid_values[k]- elevation_centroid_values[k] + elevation_edge_values[k3+1]);
    1012                   //stage_edge_values[k3+2] = min(stage_centroid_values[k], stage_centroid_values[k]- elevation_centroid_values[k] + elevation_edge_values[k3+2]);
    1013                   // Isolated wet cell -- constant stage extrapolation
    1014                   stage_edge_values[k3]   = stage_centroid_values[k];
    1015                   stage_edge_values[k3+1] = stage_centroid_values[k];
    1016                   stage_edge_values[k3+2] = stage_centroid_values[k];
    1017               }else{
    1018                   // Cell with one wet neighbour.
    1019                   // Find which neighbour is wet [if k!=k0, then k0 is wet]
    1020                   if(k!=k0){
    1021                     ktmp=k0;
    1022                     ii=0;
    1023                     xtmp = x0;
    1024                     ytmp = y0;
    1025                   }else if(k!=k1){
    1026                     ktmp=k1;
    1027                     ii=1;
    1028                     xtmp = x1;
    1029                     ytmp = y1;
    1030                   }else if(k!=k2){
    1031                     ktmp=k2;
    1032                     ii=2;
    1033                     xtmp = x2;
    1034                     ytmp = y2;
    1035                   }else{
    1036                       printf("ERROR in extrapolation routine: Should have had one ki == k \n");
    1037                       report_python_error(AT, "Negative triangle area");
    1038                       return -1;
    1039                   }
    1040 
    1041                   //if(stage_centroid_values[k]>elevation_centroid_values[k]+minimum_allowed_height){
    1042                   if(k_wetdry==0){
    1043                   //    // Cell is wet
    1044                   //    //if(count_wet_neighbours[ktmp]>0){
    1045                   //    //    stage_edge_values[k3]= stage_centroid_values[k];
    1046                   //    //    stage_edge_values[k3+1]= stage_centroid_values[k];
    1047                   //    //    stage_edge_values[k3+2]= stage_centroid_values[k];
    1048                   //    //}else{
    1049                   //        // Constant depth extrapolation
    1050                       //if(stage_centroid_values[k]<stage_centroid_values[ktmp]){
    1051                           //stage_edge_values[k3]= stage_centroid_values[k]-elevation_centroid_values[k]+elevation_edge_values[k3];
    1052                           //stage_edge_values[k3+1]= stage_centroid_values[k] -elevation_centroid_values[k]+elevation_edge_values[k3+1];
    1053                           //stage_edge_values[k3+2]= stage_centroid_values[k] -elevation_centroid_values[k]+elevation_edge_values[k3+2];
    1054                       stage_edge_values[k3] = stage_centroid_values[k];
    1055                       stage_edge_values[k3+1] = stage_centroid_values[k];
    1056                       stage_edge_values[k3+2] = stage_centroid_values[k];
    1057                       //}else{
    1058                       //    stage_edge_values[k3] = stage_centroid_values[ktmp];
    1059                       //    stage_edge_values[k3+1] = stage_centroid_values[ktmp];
    1060                       //    stage_edge_values[k3+2] = stage_centroid_values[ktmp];
    1061                       //    //stage_edge_values[k3]= max(stage_centroid_values[k]-elevation_centroid_values[k]+elevation_edge_values[k3], stage_centroid_values[ktmp]);
    1062                       //    //stage_edge_values[k3+1]= max(stage_centroid_values[k] -elevation_centroid_values[k]+elevation_edge_values[k3+1],stage_centroid_values[ktmp]);
    1063                       //    //stage_edge_values[k3+2]= max(stage_centroid_values[k] -elevation_centroid_values[k]+elevation_edge_values[k3+2], stage_centroid_values[ktmp]);
    1064                       //}
    1065 
    1066                   }else{
    1067                      
    1068                   //    // This cell is 'dry'.
    1069                   //    // Extrapolate using the neighbouring wet cell if appropriate
    1070                   //    // This needs to preserve a wet-dry interface
    1071                       stage_edge_values[k3]= stage_edge_values[3*ktmp+neighbour_edges[3*ktmp+ii]];//max(stage_edge_values[3*ktmp+neighbour_edges[3*ktmp+i]], elevation_edge_values[k3]) ;
    1072                       stage_edge_values[k3+1]= stage_edge_values[3*ktmp+neighbour_edges[3*ktmp+ii]];//max(stage_edge_values[3*ktmp+neighbour_edges[3*ktmp+i]], elevation_edge_values[k3+1]);
    1073                       stage_edge_values[k3+2]= stage_edge_values[3*ktmp+neighbour_edges[3*ktmp+ii]];//max(stage_edge_values[3*ktmp+neighbour_edges[3*ktmp +i]], elevation_edge_values[k3+2]);
    1074                   //    stage_edge_values[k3]= stage_centroid_values[ktmp] ;
    1075                   //    stage_edge_values[k3+1]= stage_centroid_values[ktmp];
    1076                   //    stage_edge_values[k3+2]= stage_centroid_values[ktmp];
    1077                   //    stage_edge_values[k3]= stage_centroid_values[k] ;
    1078                   //    stage_edge_values[k3+1]= stage_centroid_values[k];
    1079                   //    stage_edge_values[k3+2]= stage_centroid_values[k];
    1080                   }
    1081 
    1082               }
     1008            //if((cell_wetness_scale[k]==1.0) | (k==k0)&&(k==k1)&&(k==k2)){
     1009              stage_edge_values[k3]   = stage_centroid_values[k];
     1010              stage_edge_values[k3+1] = stage_centroid_values[k];
     1011              stage_edge_values[k3+2] = stage_centroid_values[k];
     1012            //}else{
     1013              // Dry cell with a single 'wet' neighbour
     1014
     1015              // Compute indices of neighbouring wet cell / edge
     1016              //ktmp = k0*(1-(k0==k)) + k1*(1-(k1==k)) + k2*(1-(k2==k));
     1017              //ii = 0*(1-(k0==k)) + 1*(1-(k1==k)) + 2*(1-(k2==k));
     1018              //if((ii<0)|(ii>2)){
     1019              //  printf("Invalid ii value \n");
     1020              //}
     1021           
     1022              //ktmp = 3*ktmp+neighbour_edges[3*k+ii];
     1023              //stage_edge_values[k3]= stage_edge_values[ktmp];//max(stage_edge_values[3*ktmp+neighbour_edges[3*ktmp+i]], elevation_edge_values[k3]) ;
     1024              //stage_edge_values[k3+1]= stage_edge_values[ktmp];//max(stage_edge_values[3*ktmp+neighbour_edges[3*ktmp+i]], elevation_edge_values[k3+1]);
     1025              //stage_edge_values[k3+2]= stage_edge_values[ktmp];//max(stage_edge_values[3*ktmp+neighbour_edges[3*ktmp +i]], elevation_edge_values[k3+2]);
     1026
     1027            //}   
     1028              //if( (k==k0) & (k==k1) & (k==k2)){
     1029              //    // Isolated wet cell  -- use constant depth extrapolation for stage
     1030              //    //stage_edge_values[k3]   = stage_centroid_values[k]- elevation_centroid_values[k] + elevation_edge_values[k3];
     1031              //    //stage_edge_values[k3+1] = stage_centroid_values[k]- elevation_centroid_values[k] + elevation_edge_values[k3+1];
     1032              //    //stage_edge_values[k3+2] = stage_centroid_values[k]- elevation_centroid_values[k] + elevation_edge_values[k3+2];
     1033
     1034              //    // Isolated wet cell -- constant depth extrapolation, but don't let water flow uphill
     1035              //    //stage_edge_values[k3]   = min(stage_centroid_values[k], stage_centroid_values[k]- elevation_centroid_values[k] + elevation_edge_values[k3]);
     1036              //    //stage_edge_values[k3+1] = min(stage_centroid_values[k], stage_centroid_values[k]- elevation_centroid_values[k] + elevation_edge_values[k3+1]);
     1037              //    //stage_edge_values[k3+2] = min(stage_centroid_values[k], stage_centroid_values[k]- elevation_centroid_values[k] + elevation_edge_values[k3+2]);
     1038              //    // Isolated wet cell -- constant stage extrapolation
     1039              //    stage_edge_values[k3]   = stage_centroid_values[k];
     1040              //    stage_edge_values[k3+1] = stage_centroid_values[k];
     1041              //    stage_edge_values[k3+2] = stage_centroid_values[k];
     1042              //}else{
     1043              //    // Cell with one wet neighbour.
     1044              //    // Find which neighbour is wet [if k!=k0, then k0 is wet]
     1045              //    if(k!=k0){
     1046              //      ktmp=k0;
     1047              //      ii=0;
     1048              //      xtmp = x0;
     1049              //      ytmp = y0;
     1050              //    }else if(k!=k1){
     1051              //      ktmp=k1;
     1052              //      ii=1;
     1053              //      xtmp = x1;
     1054              //      ytmp = y1;
     1055              //    }else if(k!=k2){
     1056              //      ktmp=k2;
     1057              //      ii=2;
     1058              //      xtmp = x2;
     1059              //      ytmp = y2;
     1060              //    }else{
     1061              //        printf("ERROR in extrapolation routine: Should have had one ki == k \n");
     1062              //        report_python_error(AT, "Negative triangle area");
     1063              //        return -1;
     1064              //    }
     1065
     1066              //    //if(stage_centroid_values[k]>elevation_centroid_values[k]+minimum_allowed_height){
     1067              //    if(k_wetdry==0){
     1068              //    //    // Cell is wet
     1069              //    //    //if(count_wet_neighbours[ktmp]>0){
     1070              //    //    //    stage_edge_values[k3]= stage_centroid_values[k];
     1071              //    //    //    stage_edge_values[k3+1]= stage_centroid_values[k];
     1072              //    //    //    stage_edge_values[k3+2]= stage_centroid_values[k];
     1073              //    //    //}else{
     1074              //    //        // Constant depth extrapolation
     1075              //        //if(stage_centroid_values[k]<stage_centroid_values[ktmp]){
     1076              //          //  stage_edge_values[k3]= stage_centroid_values[k]-elevation_centroid_values[k]+elevation_edge_values[k3];
     1077              //          //  stage_edge_values[k3+1]= stage_centroid_values[k] -elevation_centroid_values[k]+elevation_edge_values[k3+1];
     1078              //          //  stage_edge_values[k3+2]= stage_centroid_values[k] -elevation_centroid_values[k]+elevation_edge_values[k3+2];
     1079              //        stage_edge_values[k3] = stage_centroid_values[k];
     1080              //        stage_edge_values[k3+1] = stage_centroid_values[k];
     1081              //        stage_edge_values[k3+2] = stage_centroid_values[k];
     1082              //        //}else{
     1083              //           // stage_edge_values[k3] = stage_centroid_values[ktmp];
     1084              //           // stage_edge_values[k3+1] = stage_centroid_values[ktmp];
     1085              //           // stage_edge_values[k3+2] = stage_centroid_values[ktmp];
     1086              //        //    //stage_edge_values[k3]= max(stage_centroid_values[k]-elevation_centroid_values[k]+elevation_edge_values[k3], stage_centroid_values[ktmp]);
     1087              //        //    //stage_edge_values[k3+1]= max(stage_centroid_values[k] -elevation_centroid_values[k]+elevation_edge_values[k3+1],stage_centroid_values[ktmp]);
     1088              //        //    //stage_edge_values[k3+2]= max(stage_centroid_values[k] -elevation_centroid_values[k]+elevation_edge_values[k3+2], stage_centroid_values[ktmp]);
     1089              //        //}
     1090
     1091              //    }else{
     1092              //       
     1093              //    //    // This cell is 'dry'.
     1094              //    //    // Extrapolate using the neighbouring wet cell if appropriate
     1095              //    //    // This needs to preserve a wet-dry interface
     1096              //        //ktmp = 3*ktmp+neighbour_edges[3*k+ii];
     1097              //        //stage_edge_values[k3]= stage_edge_values[ktmp];//max(stage_edge_values[3*ktmp+neighbour_edges[3*ktmp+i]], elevation_edge_values[k3]) ;
     1098              //        //stage_edge_values[k3+1]= stage_edge_values[ktmp];//max(stage_edge_values[3*ktmp+neighbour_edges[3*ktmp+i]], elevation_edge_values[k3+1]);
     1099              //        //stage_edge_values[k3+2]= stage_edge_values[ktmp];//max(stage_edge_values[3*ktmp+neighbour_edges[3*ktmp +i]], elevation_edge_values[k3+2]);
     1100              //    //    stage_edge_values[k3]= stage_centroid_values[ktmp] ;
     1101              //    //    stage_edge_values[k3+1]= stage_centroid_values[ktmp];
     1102              //    //    stage_edge_values[k3+2]= stage_centroid_values[ktmp];
     1103              //        stage_edge_values[k3]= stage_centroid_values[k] ;
     1104              //        stage_edge_values[k3+1]= stage_centroid_values[k];
     1105              //        stage_edge_values[k3+2]= stage_centroid_values[k];
     1106              //    }
     1107              //
     1108              //}
    10831109              // First order momentum / velocity extrapolation
    10841110              //stage_edge_values[k3]    = stage_centroid_values[k];
    … …  
    11471173          // Limit the gradient
    11481174          // This is more complex for stage than other quantities
    1149           if((count_wet_neighbours[k]>0)){ //&&(stage_centroid_values[k] > elevation_centroid_values[k])){
     1175          if((count_wet_neighbours[k]>0)){//&&(k_wetdry==0)){
    11501176              limit_gradient(dqv, qmin, qmax, beta_tmp);
    11511177              stage_edge_values[k3+0] = stage_centroid_values[k] + dqv[0];
    … …  
    11531179              stage_edge_values[k3+2] = stage_centroid_values[k] + dqv[2];
    11541180          }else{
    1155           //if(count_wet_neighbours[k]==0){
    1156               //weight=max(neighbour_wetness_scale[k] - (stage_centroid_values[k]-elevation_centroid_values[k]), 0.);
    1157               //weight/=(max_elevation_edgevalue[k] - elevation_centroid_values[k]);
    1158               //weight=min(weight, 1.0);
    1159               //weight==0;
    1160               //weight=0.;
    1161               // 'Shallow flow on steep slope'
     1181          ////if(count_wet_neighbours[k]==0){
     1182          //    //weight=max(cell_wetness_scale[k] - (stage_centroid_values[k]-elevation_centroid_values[k]), 0.);
     1183          //    //weight/=(max_elevation_edgevalue[k] - elevation_centroid_values[k]);
     1184          //    //weight=min(weight, 1.0);
     1185          //    //weight==0;
     1186          //    //weight=0.;
     1187          //    // 'Shallow flow on steep slope'
     1188          //    //limit_gradient(dqv, qmin, qmax, beta_tmp);
    11621189              stage_edge_values[k3+0] =stage_centroid_values[k] + dqv[0];
    11631190              stage_edge_values[k3+1] =stage_centroid_values[k] + dqv[1];
    11641191              stage_edge_values[k3+2] =stage_centroid_values[k] + dqv[2];
    11651192
    1166               // There exists a neighbouring bed cell with elevation > minimum
    1167               // neighbouring stage value
    1168 
    1169               // We have to be careful in this situation. Limiting the stage gradient can
    1170               // cause problems for shallow flows down steep slopes (rainfall-runoff type problems)
    1171               // Previously 'balance_deep_and_shallow' was used to deal with this issue
    1172               //for(i=0; i<3; i++){
    1173               //    stage_edge_values[k3+i] = weight*stage_edge_values[k3+i] +
    1174               //                             (1.0-weight)*(stage_centroid_values[k] -elevation_centroid_values[k]
    1175               //                                           +elevation_edge_values[k3+i]);
    1176               //}
     1193          //    // There exists a neighbouring bed cell with elevation > minimum
     1194          //    // neighbouring stage value
     1195
     1196          //    // We have to be careful in this situation. Limiting the stage gradient can
     1197          //    // cause problems for shallow flows down steep slopes (rainfall-runoff type problems)
     1198          //    // Previously 'balance_deep_and_shallow' was used to deal with this issue
     1199          //    //for(i=0; i<3; i++){
     1200          //    //    stage_edge_values[k3+i] = weight*stage_edge_values[k3+i] +
     1201          //    //                             (1.0-weight)*(stage_centroid_values[k] -elevation_centroid_values[k]
     1202          //    //                                           +elevation_edge_values[k3+i]);
     1203          //    //}
    11771204          }
    11781205           
    … …  
    12141241          //if(stagemin>=bedmax){
    12151242          if((count_wet_neighbours[k]>0) &&
    1216             (stage_centroid_values[k]>(minimum_allowed_height+elevation_centroid_values[k]))){
     1243             (cell_wetness_scale[k]==1.0)){
     1244            //(stage_centroid_values[k]>(minimum_allowed_height+elevation_centroid_values[k]))){
    12171245            limit_gradient(dqv, qmin, qmax, beta_tmp);
    12181246          }else{
    … …  
    12661294          //if(stagemin>=bedmax){
    12671295          if((count_wet_neighbours[k]>0)&&
    1268             (stage_centroid_values[k]>(minimum_allowed_height+elevation_centroid_values[k]))){
     1296             (cell_wetness_scale[k]==1.0)){
     1297            //(stage_centroid_values[k]>(minimum_allowed_height+elevation_centroid_values[k]))){
    12691298            limit_gradient(dqv, qmin, qmax, beta_tmp);
    12701299          }else{
    … …  
    14571486  //    // Hack for 'dry' cells
    14581487  //    // Make sure edge value is not greater than neighbour edge value
    1459   //    if(stage_centroid_values[k] - elevation_centroid_values[k] < minimum_allowed_height+epsilon){
     1488  //    //if(stage_centroid_values[k] - elevation_centroid_values[k] < minimum_allowed_height+epsilon){
     1489  //    if(cell_wetness_scale[k]==0.0){
    14601490  //      for(i=0; i<3;i++){
    14611491  //          if(surrogate_neighbours[3*k+i] >= 0){
    14621492  //              ktmp=3*surrogate_neighbours[3*k+i] + neighbour_edges[3*k+i];
    1463   //          //if(neighbour_wetness_scale[surrogate_neighbours[3*k+i]]>0.){
    1464   //              stage_edge_values[3*k+i] = max(stage_edge_values[3*k+i], stage_edge_values[ktmp]);
     1493  //          //if(cell_wetness_scale[surrogate_neighbours[3*k+i]]>0.){
     1494  //              stage_edge_values[3*k+i] = min(stage_edge_values[3*k+i], stage_edge_values[ktmp]);
    14651495  //          }
    14661496  //      }
    … …  
    15181548  free(ymom_centroid_store);
    15191549  free(stage_centroid_store);
    1520   free(min_elevation_edgevalue);
     1550  //free(min_elevation_edgevalue);
    15211551  free(max_elevation_edgevalue);
    1522   free(neighbour_wetness_scale);
     1552  free(cell_wetness_scale);
    15231553  free(count_wet_neighbours);
    15241554

  • trunk/anuga_work/development/gareth/tests/runup_sinusoid/runup_sinusoid.py

    r8547 r8549  
    3737
    3838def stagefun(x,y):
    39     stge=-0.2*scale_me # +0.01*(x>0.9)
     39    stge=-0.2*scale_me #+0.01*(x>0.9)
    4040    #topo=topography(x,y)
    4141    return stge#*(stge>topo) + (topo)*(stge<=topo)
    … …  
    7272# Associate boundary tags with boundary objects
    7373#----------------------------------------------
    74 domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom':Br})
     74domain.set_boundary({'left': Br, 'right': Bd, 'top': Br, 'bottom':Br})
    7575
    7676#------------------------------
    … …  
    7878#------------------------------
    7979
    80 for t in domain.evolve(yieldstep=0.2,finaltime=20.0):
     80for t in domain.evolve(yieldstep=0.2,finaltime=40.0):
    8181    print domain.timestepping_statistics()
    8282    xx = domain.quantities['xmomentum'].centroid_values

  • trunk/anuga_work/development/gareth/tests/shallow_steep_slope/channel_SU_sparse.py

    r8547 r8549  
    3030#------------------------------------------------------------------------------
    3131points, vertices, boundary = rectangular_cross(10, 10,
    32 len1=100.0, len2=100.0) # Mesh
     32                                len1=100.0, len2=100.0) # Mesh
    3333domain = Domain(points, vertices, boundary) # Create domain
    3434domain.set_name('channel_SU_2_v2') # Output name
    35 domain.set_store_vertices_uniquely(True)
    36 domain.CFL=1.0
     35#domain.set_store_vertices_uniquely(True)
     36#domain.CFL=1.0
    3737#------------------------------------------------------------------------------
    3838# Setup initial conditions
    3939#------------------------------------------------------------------------------
    4040def topography(x, y):
    41         return -x/10. + numpy.sin(x/10.) +abs(y-50.)/100. # linear bed slope
     41        return -x/10. #+ numpy.sin(x/10.) +abs(y-50.)/100. # linear bed slope
    4242
    4343def stagetopo(x,y):
    … …  
    5353domain.set_quantity('friction', 0.03) # Constant friction
    5454domain.set_quantity('stage', stagetopo)
     55
    5556#------------------------------------------------------------------------------
    5657# Setup boundary conditions
    … …  
    6970domain.set_boundary({'left': Br, 'right': Bo, 'top': Br, 'bottom': Br})
    7071
     72
     73myindex=((domain.centroid_coordinates[:,0] - 50.)**2 + (domain.centroid_coordinates[:,1] - 50.)**2).argmin()
     74
    7175#------------------------------------------------------------------------------
    7276# Evolve system through time
    7377#------------------------------------------------------------------------------
    7478for t in domain.evolve(yieldstep=4.0, finaltime=400.0):
    75         print domain.timestepping_statistics()
     79    print domain.timestepping_statistics()
     80    #print domain.quantities['stage'].centroid_values[myindex] - domain.quantities['elevation'].centroid_values[myindex]
     81    s3 = domain.get_flow_through_cross_section([[30., 0.0], [30., 100.]])
     82    s4 = domain.get_flow_through_cross_section([[32., 0.0], [32., 100.]])
     83    s5 = domain.get_flow_through_cross_section([[34., 0.0], [34., 100.]])
     84    s2 = domain.get_flow_through_cross_section([[45., 0.0], [45., 100.]])
     85    s1 = domain.get_flow_through_cross_section([[53., 0.0], [53., 100.]])
     86    s0 = domain.get_flow_through_cross_section([[60., 0.0], [60., 100.]])
     87    print 'Xsectional flow:', s0, s1, s2, s3, s4, s5
     88   
Note: See TracChangeset for help on using the changeset viewer.