Changeset 8884


Ignore:
Timestamp:
Jun 7, 2013, 3:49:59 PM (11 years ago)
Author:
davies
Message:

Updates to bal_dev

Location:
trunk/anuga_work/development/gareth
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain.py

    r8866 r8884  
    6666        self.set_CFL(1.0)
    6767        self.set_use_kinematic_viscosity(False)
    68         self.timestepping_method='rk2'#'rk2'#'rk2'#'euler'#'rk2'
     68        self.timestepping_method='rk2'#'euler'#'rk3'#'euler'#'rk2'
    6969        # The following allows storage of the negative depths associated with this method
    7070        self.minimum_storable_height=-99999999999.0
    … …  
    7474        self.extrapolate_velocity_second_order=True
    7575
    76         self.beta_w=0.0
     76        self.beta_w=1.0
    7777        self.beta_w_dry=0.0
    78         self.beta_uh=0.0
     78        self.beta_uh=1.0
    7979        self.beta_uh_dry=0.0
    80         self.beta_vh=0.0
     80        self.beta_vh=1.0
    8181        self.beta_vh_dry=0.0
    8282

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

    r8867 r8884  
    234234  // Flux formulas
    235235  flux_left[0] = u_left*h_left;
    236   //if(hc>h0){
     236  //if(hc>h0 & hc_n > h0){
     237  //if(w_left>z+h0 & w_right>z+h0){
    237238    flux_left[1] = u_left*uh_left; //+ 0.5*g*h_left*h_left;
    238239    flux_left[2] = u_left*vh_left;
    … …  
    243244
    244245  flux_right[0] = u_right*h_right;
    245   //if(hc_n>h0){
     246  //if(w_left>z+h0 & w_right>z+h0){
    246247    flux_right[1] = u_right*uh_right ; //+ 0.5*g*h_right*h_right;
    247248    flux_right[2] = u_right*vh_right;
    … …  
    253254  // Flux computation
    254255  denom = s_max - s_min;
    255   //if (denom < -epsilon)
     256  //if (denom < epsilon)
    256257  if (0==1)
    257258  {
    … …  
    275276
    276277      // Standard smoothing term
    277       //edgeflux[i] += (s_max*s_min)*(q_right_rotated[i] - q_left_rotated[i]);
     278      edgeflux[i] += (s_max*s_min)*(q_right_rotated[i] - q_left_rotated[i]);
    278279
    279280      // Standard smoothing term, scaled
    … …  
    282283      // Adjustment to the scheme by Kurganov and Lin 2007 Communications in Computational
    283284      // Physics 2:141-163
    284       uint = (s_max*q_right_rotated[i] - s_min*q_left_rotated[i] - (flux_right[i] - flux_left[i]))*inverse_denominator;
    285       t1 = (q_right_rotated[i] - uint);
    286       t2 = (-q_left_rotated[i] + uint);
    287       t3 = _minmod(t1,t2);
    288       edgeflux[i] += (s_max*s_min)*(q_right_rotated[i] - q_left_rotated[i]-t3);
     285      //uint = (s_max*q_right_rotated[i] - s_min*q_left_rotated[i] - (flux_right[i] - flux_left[i]))*inverse_denominator;
     286      //t1 = (q_right_rotated[i] - uint);
     287      //t2 = (-q_left_rotated[i] + uint);
     288      //t3 = _minmod(t1,t2);
     289      //edgeflux[i] += (s_max*s_min)*(q_right_rotated[i] - q_left_rotated[i]-t3);
    289290
    290291      // test this
    … …  
    457458            }
    458459           
    459             // GD HACK -- although I think this is a common technique
    460             //If one centroid is dry, then extrapolate its edge values from the neighbouring centroid,
    461             // unless the local centroid value is smaller
    462             //if(n>=0){
    463             //    if((hc<=H0)&&(hc_n>H0)){
    464             //        ql[0]=max(min(qr[0],stage_centroid_values[k]),zl);
    465             //    }
    466             //    if((hc_n<=H0)&&(hc>H0)){
    467             //        qr[0]=max(min(ql[0],stage_centroid_values[n]),zr);
    468             //    }
    469 
    470             //}else{
    471             //    // Treat the boundary case
    472             //    // ??
    473             //}
    474460
    475461            smooth=1.0;
    … …  
    500486            // bedslope_work contains all gravity related terms -- weighting of
    501487            // conservative and non-conservative versions.
     488
     489            if(ql[0]>zl){
     490               tmp=1.0;
     491            }else{
     492               tmp=1.;
     493            }
    502494            if(hc> h0){
    503495            //if(stage_centroid_values[k] > count_wet_edges[k]){
    504               bedslope_work = length*(g*(hc*ql[0]-0.5*max(ql[0]-zl,0.)*(ql[0]-zl)) + pressure_flux);
     496              bedslope_work = length*(g*(hc*ql[0]-0.5*max(ql[0]-zl,0.)*(ql[0]-zl))*tmp + pressure_flux);
    505497              //bedslope_work = length*(g*(hc*stage_centroid_values[k]*0.-
    506498              //                           0.5*max(stage_centroid_values[k]-zl,0.)*(stage_centroid_values[k]-zl))
    … …  
    509501              pressuregrad_store[ki] = bedslope_work;
    510502            }else{
    511               // NOTE: When hc<h0, pressure_flux is entirely computed from the
    512               // other side of the edge. Hence, we can't necessarily satisfy
    513               // well-balancing by just assuming that the 2nd and 3rd terms in
    514               // bedslope_work cancel. So, we force it here -- equivalent to
    515               // using the water surface gravity term
    516               bedslope_work = length*(g*(hc*ql[0]-0.0*max(ql[0]-zl,0.)*(ql[0]-zl)) + 0.0*pressure_flux);
     503              // NOTE: When hc<h0, local contribution to the pressure_flux is
     504              // entirely computed from the other side of the edge. Hence, we
     505              // can't necessarily satisfy well-balancing by just assuming that
     506              // the 2nd and 3rd terms in bedslope_work cancel. So, we force it
     507              // here -- equivalent to using the water surface gravity term
     508              bedslope_work = length*(g*(hc*ql[0]*tmp-0.0*max(ql[0]-zl,0.)*(ql[0]-zl))*tmp + 0.0*pressure_flux);
    517509              //bedslope_work = length*(g*(hc*stage_centroid_values[k]-0.0*max(ql[0]-zl,0.)*(ql[0]-zl)) + 0.0*pressure_flux);
    518510              //bedslope_work+= (1.0-tmp)*length*g*hc*ql[0]; // Non-conservative water surface slope
    … …  
    532524                edgeflux_store[nm3 + 2 ] = edgeflux[2];
    533525
     526                if(qr[0]>zr){
     527                   tmp=1.0;
     528                }else{
     529                   tmp=1.;
     530                }
    534531                if(hc_n> h0){
    535532                //if(stage_centroid_values[n] > count_wet_edges[n]){
    536                   bedslope_work = length*(g*(hc_n*qr[0]-0.5*max(qr[0]-zr,0.)*(qr[0]-zr)) + pressure_flux);
     533                  bedslope_work = length*(g*(hc_n*qr[0]*tmp-0.5*max(qr[0]-zr,0.)*(qr[0]-zr))*tmp + pressure_flux);
    537534                  //bedslope_work = length*(g*(hc_n*stage_centroid_values[n]*0.-0.5*max(stage_centroid_values[n]-zr,0.)*
    538535                  //                                (stage_centroid_values[n]-zr)) + pressure_flux);
    … …  
    545542                  // bedslope_work cancel. So, we force it here -- equivalent to
    546543                  // using the water surface gravity term
    547                   bedslope_work = length*(g*(hc_n*qr[0]-0.0*max(qr[0]-zr,0.)*(qr[0]-zr)) + 0.0*pressure_flux);
     544                  bedslope_work = length*(g*(hc_n*qr[0]*tmp-0.0*max(qr[0]-zr,0.)*(qr[0]-zr))*tmp + 0.0*pressure_flux);
    548545                  //bedslope_work = length*(g*(hc_n*stage_centroid_values[n]-0.0*max(qr[0]-zr,0.)*(qr[0]-zr)) + 0.0*pressure_flux);
    549546                  //bedslope_work+= (1.-tmp)*length*g*hc_n*qr[0]; 
    … …  
    778775  // distance between the bed_centroid_value and the max bed_edge_value of
    779776  // every triangle.
    780   double minimum_relative_height=0.1;
     777  double minimum_relative_height=0.05;
    781778  int mass_added=0;
    782779
    … …  
    786783      hc = wc[k] - zc[k];
    787784      // Definine the maximum bed edge value on triangle k.
    788       bmax = 0.5*max((zv[3*k]+zv[3*k+1]),max((zv[3*k+1]+zv[3*k+2]),(zv[3*k+2]+zv[3*k])));
     785      //bmax = 0.5*max((zv[3*k]+zv[3*k+1]),max((zv[3*k+1]+zv[3*k+2]),(zv[3*k+2]+zv[3*k])));
     786      bmax = max(zv[3*k],max(zv[3*k+2],zv[3*k+1]));
    789787
    790788      //if (hc < max(minimum_relative_height*(bmax-zc[k]), minimum_allowed_height)) {
    … …  
    796794      //if (hc < minimum_allowed_height*2.0 | hc < minimum_relative_height*(bmax-zc[k]) ){
    797795      if (hc < minimum_allowed_height*2.0 ){
     796      //if (hc < minimum_allowed_height*2.0 ){
    798797            // Set momentum to zero and ensure h is non negative
    799             xmomc[k] = 0.;//xmomc[k]*hc*max(hc-minimum_allowed_height,0.)/(2.0*minimum_allowed_height*minimum_allowed_height);// 0.0;
    800             ymomc[k] = 0.;//ymomc[k]*hc*max(hc-minimum_allowed_height,0.)/(2.0*minimum_allowed_height*minimum_allowed_height);// 0.0;
     798            //xmomc[k] = 0.;//xmomc[k]*hc*max(hc-minimum_allowed_height,0.)/(2.0*minimum_allowed_height*minimum_allowed_height);// 0.0;
     799            //ymomc[k] = 0.;//ymomc[k]*hc*max(hc-minimum_allowed_height,0.)/(2.0*minimum_allowed_height*minimum_allowed_height);// 0.0;
    801800            //xmomc[k] = xmomc[k]*hc*max(hc-minimum_allowed_height,0.)/(2.0*minimum_allowed_height*minimum_allowed_height);// 0.0;
    802801            //ymomc[k] = ymomc[k]*hc*max(hc-minimum_allowed_height,0.)/(2.0*minimum_allowed_height*minimum_allowed_height);// 0.0;
    … …  
    809808             //bmin =0.5*bmin + 0.5*min(zv[3*k],min(zv[3*k+1],zv[3*k+2]));
    810809             //bmin =min(zv[3*k],min(zv[3*k+1],zv[3*k+2])) -minimum_allowed_height;
    811              bmin=zc[k];//-1.0e-03;
     810             bmin=zc[k]-minimum_allowed_height;//-1.0e-03;
    812811             // Minimum allowed stage = bmin
    813812
    … …  
    10061005  }
    10071006
     1007
     1008  // Alternative 'PROTECT' step
     1009  for(k=0; k<number_of_elements;k++){
     1010    if(count_wet_neighbours[k]<=3){
     1011        bedmax = max(elevation_centroid_values[k],
     1012                     max(elevation_vertex_values[3*k],
     1013                         max(elevation_vertex_values[3*k+1],
     1014                             elevation_vertex_values[3*k+2])));
     1015        if(cell_wetness_scale[k]==0. ||
     1016            stage_centroid_values[k] - elevation_centroid_values[k] < 0.05*(bedmax-elevation_centroid_values[k])){
     1017            xmom_centroid_store[k]=0.;
     1018            xmom_centroid_values[k]=0.;
     1019            ymom_centroid_store[k]=0.;
     1020            ymom_centroid_values[k]=0.;
     1021        }
     1022    }
     1023  }
    10081024 
    10091025  // First treat all 'fully wet' cells, then treat 'partially dry' cells
    … …  
    11281144                           max(elevation_centroid_values[k1],
    11291145                               elevation_centroid_values[k2])));
    1130           bedmin = min(elevation_centroid_values[k],
    1131                        min(elevation_centroid_values[k0],
    1132                            min(elevation_centroid_values[k1],
    1133                                elevation_centroid_values[k2])));
    1134           //stagemin = min(max(stage_centroid_values[k], elevation_centroid_values[k]),
    1135           //               min(max(stage_centroid_values[k0], elevation_centroid_values[k0]),
    1136           //                   min(max(stage_centroid_values[k1], elevation_centroid_values[k1]),
    1137           //                       max(stage_centroid_values[k2], elevation_centroid_values[k2]))));
     1146          //bedmin = min(elevation_centroid_values[k],
     1147          //             min(elevation_centroid_values[k0],
     1148          //                 min(elevation_centroid_values[k1],
     1149          //                     elevation_centroid_values[k2])));
     1150          stagemin = min(max(stage_centroid_values[k], elevation_centroid_values[k]),
     1151                         min(max(stage_centroid_values[k0], elevation_centroid_values[k0]),
     1152                             min(max(stage_centroid_values[k1], elevation_centroid_values[k1]),
     1153                                 max(stage_centroid_values[k2], elevation_centroid_values[k2]))));
    11381154
    11391155          // Get the auxiliary triangle's vertex coordinates
    … …  
    11631179           
    11641180          //// Treat triangles with zero or 1 wet neighbours (area2 <=0.)
    1165           if ((area2 <= 0.))// || (cell_wetness_scale[k]==0.)) //|(count_wet_neighbours[k]==0))
     1181          if ((area2 <= 0.) )// || (cell_wetness_scale[k]==0.)) //|(count_wet_neighbours[k]==0))
    11661182          {
    11671183
    … …  
    12171233
    12181234            //}   
    1219               //if( (k==k0) & (k==k1) & (k==k2)){
     1235              if( (k==k0) & (k==k1) & (k==k2)){
    12201236              //    // Isolated wet cell  -- use constant depth extrapolation for stage
    12211237              //    stage_edge_values[k3]   = stage_centroid_values[k]- elevation_centroid_values[k] + elevation_edge_values[k3];
    … …  
    12281244              //    //stage_edge_values[k3+2] = min(stage_centroid_values[k], stage_centroid_values[k]- elevation_centroid_values[k] + elevation_edge_values[k3+2]);
    12291245              //    // Isolated wet cell -- constant stage extrapolation
    1230               //    stage_edge_values[k3]   = stage_centroid_values[k];
    1231               //    stage_edge_values[k3+1] = stage_centroid_values[k];
    1232               //    stage_edge_values[k3+2] = stage_centroid_values[k];
    1233               //}else{
     1246                  stage_edge_values[k3]   = stage_centroid_values[k];
     1247                  stage_edge_values[k3+1] = stage_centroid_values[k];
     1248                  stage_edge_values[k3+2] = stage_centroid_values[k];
     1249              }else{
    12341250              //    // Cell with one wet neighbour.
    12351251              //    // Find which neighbour is wet [if k!=k0, then k0 is wet]
    1236               //    if(k!=k0){
    1237               //      ktmp=k0;
    1238               //      ii=0;
    1239               //      xtmp = x0;
    1240               //      ytmp = y0;
    1241               //    }else if(k!=k1){
    1242               //      ktmp=k1;
    1243               //      ii=1;
    1244               //      xtmp = x1;
    1245               //      ytmp = y1;
    1246               //    }else if(k!=k2){
    1247               //      ktmp=k2;
    1248               //      ii=2;
    1249               //      xtmp = x2;
    1250               //      ytmp = y2;
    1251               //    }else{
    1252               //        printf("ERROR in extrapolation routine: Should have had one ki == k \n");
    1253               //        report_python_error(AT, "Negative triangle area");
    1254               //        return -1;
    1255               //    }
     1252                  if(k!=k0){
     1253                    ktmp=k0;
     1254                    ii=0;
     1255                    xtmp = x0;
     1256                    ytmp = y0;
     1257                  }else if(k!=k1){
     1258                    ktmp=k1;
     1259                    ii=1;
     1260                    xtmp = x1;
     1261                    ytmp = y1;
     1262                  }else if(k!=k2){
     1263                    ktmp=k2;
     1264                    ii=2;
     1265                    xtmp = x2;
     1266                    ytmp = y2;
     1267                  }else{
     1268                      printf("ERROR in extrapolation routine: Should have had one ki == k \n");
     1269                      report_python_error(AT, "Negative triangle area");
     1270                      return -1;
     1271                  }
    12561272
    12571273              //    //if(stage_centroid_values[k]>elevation_centroid_values[k]+minimum_allowed_height){
    … …  
    12751291              //           // stage_edge_values[k3+1] = stage_centroid_values[ktmp];
    12761292              //           // stage_edge_values[k3+2] = stage_centroid_values[ktmp];
    1277               //        //    //stage_edge_values[k3]= max(stage_centroid_values[k]-elevation_centroid_values[k]+elevation_edge_values[k3], stage_centroid_values[ktmp]);
    1278               //        //    //stage_edge_values[k3+1]= max(stage_centroid_values[k] -elevation_centroid_values[k]+elevation_edge_values[k3+1],stage_centroid_values[ktmp]);
    1279               //        //    //stage_edge_values[k3+2]= max(stage_centroid_values[k] -elevation_centroid_values[k]+elevation_edge_values[k3+2], stage_centroid_values[ktmp]);
     1293                     stage_edge_values[k3]= max(stage_centroid_values[k]-elevation_centroid_values[k]+elevation_edge_values[k3], stage_centroid_values[ktmp]);
     1294                     stage_edge_values[k3+1]= max(stage_centroid_values[k] -elevation_centroid_values[k]+elevation_edge_values[k3+1],stage_centroid_values[ktmp]);
     1295                     stage_edge_values[k3+2]= max(stage_centroid_values[k] -elevation_centroid_values[k]+elevation_edge_values[k3+2], stage_centroid_values[ktmp]);
    12801296              //        //}
    12811297
    … …  
    12971313              //    }
    12981314              //
    1299               //}
     1315              }
    13001316              // First order momentum / velocity extrapolation
    1301               stage_edge_values[k3]    = stage_centroid_values[k];
    1302               stage_edge_values[k3+1]  = stage_centroid_values[k];
    1303               stage_edge_values[k3+2]  = stage_centroid_values[k];
     1317             //if(stagemin>=bedmax){
     1318             //   stage_edge_values[k3]    = stage_centroid_values[k];
     1319             //   stage_edge_values[k3+1]  = stage_centroid_values[k];
     1320             //   stage_edge_values[k3+2]  = stage_centroid_values[k];
     1321             //}else{
     1322             //   stage_edge_values[k3]= stage_centroid_values[k]-elevation_centroid_values[k]+elevation_edge_values[k3];
     1323             //   stage_edge_values[k3+1]= stage_centroid_values[k] -elevation_centroid_values[k]+elevation_edge_values[k3+1];
     1324             //   stage_edge_values[k3+2]= stage_centroid_values[k] -elevation_centroid_values[k]+elevation_edge_values[k3+2];
     1325             //}
    13041326              xmom_edge_values[k3]    = xmom_centroid_values[k];
    13051327              xmom_edge_values[k3+1]  = xmom_centroid_values[k];
    … …  
    13271349            //hfactor=hmin/(hmin + 0.004);
    13281350          //}
    1329           hfactor= min(2.0*max(hmin,0.0)/max(hc,max(0.5*(bedmax-bedmin),minimum_allowed_height*2.)), 1.0);
    13301351          //hfactor= min(2.0*max(hmin,0.0)/max(hc,max(0.5*(bedmax-bedmin),minimum_allowed_height)), 1.0);
     1352          hfactor= min(2.0*max(hmin,0.0)/max(hc,minimum_allowed_height*0.1), 1.0);
    13311353         
    13321354          //-----------------------------------
    … …  
    13921414          //if( (area2>0) ){
    13931415          //if(count_wet_neighbours[k]>0){
    1394           limit_gradient(dqv, qmin, qmax, beta_tmp);
     1416          //if(bedmax<=stagemin){
     1417             limit_gradient(dqv, qmin, qmax, beta_tmp);
     1418          //}else{
     1419          //  dqv[0]=-elevation_centroid_values[k] + elevation_edge_values[k3+0];
     1420          //  dqv[1]=-elevation_centroid_values[k] + elevation_edge_values[k3+1];
     1421          //  dqv[2]=-elevation_centroid_values[k] + elevation_edge_values[k3+2];
     1422          //}
    13951423          //printf("%lf \n ", beta_tmp);
    13961424          //}
    … …  
    14881516            // (cell_wetness_scale[k]==1.0)){
    14891517          //if(stage_centroid_values[k]>(minimum_allowed_height+max_elevation_edgevalue[k])){
     1518         
    14901519            limit_gradient(dqv, qmin, qmax, beta_tmp);
    14911520          //}else{

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

    r8880 r8884  
    1414#from swb2_domain import *
    1515#from balanced_basic import *
    16 #from balanced_dev import *
     16from balanced_dev import *
    1717#from anuga_tsunami import *
    1818
    … …  
    2222points, vertices, boundary = anuga.rectangular_cross(20,20, len1=1., len2=1.)
    2323
    24 domain=anuga.Domain(points,vertices,boundary)    # Create Domain
     24domain=Domain(points,vertices,boundary)    # Create Domain
    2525domain.set_name('runup_sinusoid_v2')                         # Output to file runup.sww
    2626#domain.set_timestepping_method('euler')
    27 domain.set_flow_algorithm('tsunami')
     27#domain.set_flow_algorithm('tsunami')
    2828#------------------
    2929# Define topography
    … …  
    3232#### Pathological
    3333#scale_me=100.0
    34 #boundary_ws=-0.2
     34#boundary_ws=-0.2#1999
    3535#init_ws=-0.2
    3636#bumpiness=50. # Higher = shorter wavelength oscillations in topography
    37 #tstep=0.002
    38 #lasttime=1.1
     37#tstep=0.2
     38#lasttime=20.1
    3939
    4040### Sensible
    … …  
    6060domain.set_quantity('friction',0.00)             # Constant friction
    6161
     62print '#-#', domain.minimum_allowed_height
    6263#def frict_change(x,y):
    6364#       return 0.2*(x>0.5)+0.1*(x<=0.5)

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

    r8865 r8884  
    4343    return stg#*(stg>topo) + topo*(stg<=topo)
    4444
    45 line1=[ [10.,0.], [10., 100.] ]
    46 Qin=20.
     45line1=[ [9.,0.], [9., 100.] ]
     46Qin=2.
    4747Inlet_operator(domain, line1,Qin)
    4848
Note: See TracChangeset for help on using the changeset viewer.