Changeset 8359


Ignore:
Timestamp:
Mar 13, 2012, 3:46:12 PM (12 years ago)
Author:
davies
Message:

balanced_dev: important change to 'protect'

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

Legend:

Unmodified
Added
Removed
  • trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain_ext.c

    r8358 r8359  
    291291    // Local variables
    292292    double max_speed, length, inv_area, zl, zr;
    293     double h0 = H0*H0; //H0*H0; // This ensures a good balance when h approaches H0.
     293    double h0 = H0*H0; // This ensures a good balance when h approaches H0.
    294294
    295295    double limiting_threshold = 10 * H0; // Avoid applying limiter below this
     
    580580
    581581  int k;
    582   double hc, wmin;
     582  double hc, bmin, bmax;
    583583  double u, v, reduced_speed;
     584
     585  // This acts like minimum_allowed height, but scales with the vertical
     586  // distance between the bed_centroid_value and the max bed_edge_value of every triangle.
     587  double minimum_relative_height=0.05;
    584588
    585589
     
    588592    for (k=0; k<N; k++) {
    589593      hc = wc[k] - zc[k];
    590 
    591       if (hc < minimum_allowed_height) {
     594      // Definine the maximum bed edge value on triangle k.
     595      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])));
     596
     597      if (hc < max(minimum_relative_height*(bmax-zc[k]), minimum_allowed_height)) {
    592598       
    593599        // Set momentum to zero and ensure h is non negative
     600        // NOTE: THIS IS IMPORTANT -- WE ARE SETTING MOMENTUM TO ZERO
    594601        xmomc[k] = 0.0;
    595602        ymomc[k] = 0.0;
     603
    596604        if (hc <= 0.0){
    597              // Definine the minimum possible stage value as the minimum bed edge value on triangle k.
    598              wmin = 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])));
    599                  wc[k] = max(wc[k], wmin) ; //max(wc[k], wmin);
    600                  wv[3*k] = max(wv[3*k], wmin);
    601                  wv[3*k+1] = max(wv[3*k+1], wmin);
    602                  wv[3*k+2] = max(wv[3*k+2], wmin);
    603              //if(wc[k]< wmin){
    604              //    //printf("Limiting wc[ %d ] = %f to %f, not %f \n", k, wc[k], max(wc[k], wmin), min(wc[k], wmin));
    605              //    wc[k] = wmin ; //max(wc[k], wmin);
    606              //    wv[3*k] = max(wv[3*k], wmin);
    607              //    wv[3*k+1] = max(wv[3*k+1], wmin);
    608              //    wv[3*k+2] = max(wv[3*k+2], wmin);
    609              //}
     605             // Definine the minimum bed edge value on triangle k.
     606             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])));
     607             // Minimum allowed stage = bmin
     608             wc[k] = max(wc[k], bmin) ;
     609             wv[3*k] = max(wv[3*k], bmin);
     610             wv[3*k+1] = max(wv[3*k+1], bmin);
     611             wv[3*k+2] = max(wv[3*k+2], bmin);
    610612        }
    611613      }
     
    780782  // Local variables
    781783  double a, b; // Gradient vector used to calculate edge values from centroids
    782   int k, k0, k1, k2, k3, k6, coord_index, i;
     784  int k, k0, k1, k2, k3, k6, coord_index, i, ii;
    783785  double x, y, x0, y0, x1, y1, x2, y2, xv0, yv0, xv1, yv1, xv2, yv2; // Vertices of the auxiliary triangle
    784786  double dx1, dx2, dy1, dy2, dxv0, dxv1, dxv2, dyv0, dyv1, dyv2, dq0, dq1, dq2, area2, inv_area2;
    785787  double dqv[3], qmin, qmax, hmin, hmax, bedmax, stagemin;
    786788  double hc, h0, h1, h2, beta_tmp, hfactor;
    787   double xmom_centroid_store[number_of_elements], ymom_centroid_store[number_of_elements], dk, dv0, dv1, dv2;
     789  double xmom_centroid_store[number_of_elements], ymom_centroid_store[number_of_elements], dk, dv0, dv1, dv2, de[3];
     790  double stage_centroid_store[number_of_elements];
    788791   
    789792  if(extrapolate_velocity_second_order==1){
     
    798801          ymom_centroid_store[k] = ymom_centroid_values[k];
    799802          ymom_centroid_values[k] = ymom_centroid_values[k]/dk;
     803
     804          // Experimental -- replace stage with a weighting of 'stage' and
     805          // 'depth', where the weights depend on the froude number **2
     806          // stage_centroid_store[k]=stage_centroid_values[k];
     807          // hfactor =(xmom_centroid_values[k]*xmom_centroid_values[k] +
     808          //          ymom_centroid_values[k]*ymom_centroid_values[k])/20.0;
     809          //          //(9.8*max(stage_centroid_values[k] - elevation_centroid_values[k],minimum_allowed_height));
     810          //stage_centroid_values[k] -= min(1.0, hfactor)*elevation_centroid_values[k];
    800811      }
    801812  }
     
    13251336      for (k=0; k<number_of_elements; k++){
    13261337          k3=3*k;
    1327           //dv0 = max(stage_edge_values[k3]-elevation_edge_values[k3],0.);
    1328           //dv1 = max(stage_edge_values[k3+1]-elevation_edge_values[k3+1],0.);
    1329           //dv2 = max(stage_edge_values[k3+2]-elevation_edge_values[k3+2],0.);
    1330 
    1331           // NOTE: dv0 etc may be negative -- we don't prevent this though,
    1332           // because if we set them to zero, then the relation between the
    1333           // centroid and edge momenta is no longer linear.
    1334           dv0 = stage_edge_values[k3]-elevation_edge_values[k3];
    1335           dv1 = stage_edge_values[k3+1]-elevation_edge_values[k3+1];
    1336           dv2 = stage_edge_values[k3+2]-elevation_edge_values[k3+2];
    1337 
    1338           //Correct centroid and edge values
     1338
     1339          // Convert energy back to stage
     1340          //stage_centroid_values[k] = stage_centroid_store[k];
     1341
     1342          //Convert velocity back to momenta
    13391343          xmom_centroid_values[k] = xmom_centroid_store[k];
    1340           xmom_edge_values[k3]=xmom_edge_values[k3]*dv0;
    1341           xmom_edge_values[k3+1]=xmom_edge_values[k3+1]*dv1;
    1342           xmom_edge_values[k3+2]=xmom_edge_values[k3+2]*dv2;
    1343          
    13441344          ymom_centroid_values[k] = ymom_centroid_store[k];
    1345           ymom_edge_values[k3]=ymom_edge_values[k3]*dv0;
    1346           ymom_edge_values[k3+1]=ymom_edge_values[k3+1]*dv1;
    1347           ymom_edge_values[k3+2]=ymom_edge_values[k3+2]*dv2;
     1345
     1346          for (i=0; i<3; i++){
     1347              //stage_edge_values[k3+i] -=(xmom_edge_values[k3+i]*xmom_edge_values[k3+i]
     1348              //                           + ymom_edge_values[k3+i]*ymom_edge_values[k3+i])*0.0;
     1349              //stage_edge_values[k3+i] += elevation_edge_values[k3+i];
     1350              //de[i] = max(stage_centroid_values[k]-elevation_centroid_values[k],minimum_allowed_height);
     1351              //for(ii=0; ii<1; ii++){
     1352              //hfactor =(xmom_edge_values[k3+i]*xmom_edge_values[k3+i] +
     1353              //          ymom_edge_values[k3+i]*ymom_edge_values[k3+i])/20.0; // 20 ~= 2*g
     1354                        //(9.8*de[i]);
     1355              //hfactor = 0.0; //max(hfactor/20.0 - 0.2, 0.0);
     1356              //stage_edge_values[k3+i] = stage_edge_values[k3+i];//+ min(hfactor, 1.0)*elevation_edge_values[k3+i];
     1357              //de[i] = max(stage_edge_values[k3+i]-elevation_edge_values[k3+i],1.0e-03);
     1358              //}
     1359
     1360              de[i] = max(stage_edge_values[k3+i]-elevation_edge_values[k3+i],0.0);
     1361              xmom_edge_values[k3+i]=xmom_edge_values[k3+i]*de[i];
     1362              ymom_edge_values[k3+i]=ymom_edge_values[k3+i]*de[i];
     1363          }
     1364
    13481365         
    13491366      }
  • trunk/anuga_work/development/gareth/tests/channel_floodplain/plotme.py

    r8355 r8359  
    44
    55# Time-index to plot outputs from
    6 index=1100
     6index=75
    77
    88#p2 = util.get_output('channel_floodplain1_bal_dev.sww', minimum_allowed_height=0.01)
  • trunk/anuga_work/development/gareth/tests/dam_break/plotme.py

    r8355 r8359  
    99p2_st=util.get_centroids(p_st)
    1010
    11 p_dev = util.get_output('dam_break_20120308_153245/dam_break.sww', 0.001)
     11p_dev = util.get_output('dam_break_20120313_153101/dam_break.sww', 0.001)
    1212p2_dev=util.get_centroids(p_dev, velocity_extrapolation=True)
    1313
  • trunk/anuga_work/development/gareth/tests/runup_sinusoid/runup_sinusoid.py

    r8358 r8359  
    1818#Setup computational domain
    1919#---------
    20 points, vertices, boundary = anuga.rectangular_cross(40,40)
     20points, vertices, boundary = anuga.rectangular_cross(40,40, len1=1., len2=1.)
    2121
    2222domain=Domain(points,vertices,boundary)    # Create Domain
     
    2525domain.set_quantities_to_be_stored({'stage': 2, 'xmomentum': 2, 'ymomentum': 2, 'elevation': 1})
    2626#domain.set_store_vertices_uniquely(True)
     27domain.minimum_allowed_height=0.001
    2728
    2829#------------------
     
    3031#------------------
    3132scale_me=1.0
     33
     34#domain.minimum_allowed_height=domain.minimum_allowed_height*scale_me # Seems needed to make the algorithms behave
     35
    3236def topography(x,y):
    3337        return (-x/2.0 +0.05*numpy.sin((x+y)*50.0))*scale_me
     
    6771# Associate boundary tags with boundary objects
    6872#----------------------------------------------
    69 domain.set_boundary({'left': Br, 'right': Bd, 'top': Br, 'bottom':Br})
     73domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom':Br})
    7074
    7175#------------------------------
     
    7579for t in domain.evolve(yieldstep=0.1,finaltime=20.0):
    7680    print domain.timestepping_statistics()
     81    xx = domain.quantities['xmomentum'].centroid_values
     82    yy = domain.quantities['ymomentum'].centroid_values
     83    dd = domain.quantities['stage'].centroid_values - domain.quantities['elevation'].centroid_values
     84    dd = (dd)*(dd>1.0e-03)+1.0e-06
     85    vv = ( (xx/dd)**2 + (yy/dd)**2 )**0.5
     86    vv = vv*(dd>1.0e-03)
     87    print 'Peak velocity is: ', vv.max(), vv.argmax()
    7788
    7889print 'Finished'
  • trunk/anuga_work/development/gareth/tests/shallow_steep_slope/channel_SU_sparse.py

    r8354 r8359  
    5757
    5858Br = anuga.Reflective_boundary(domain) # Solid reflective wall
    59 Bo= anuga.Dirichlet_boundary([-9.9369037,0.0,0]) # Outflow for steady uniform flow with a depth integrated velocity of 0.1
     59#Bo= anuga.Dirichlet_boundary([-9.9369037,0.0,0]) # Outflow for steady uniform flow with a depth integrated velocity of 0.1
     60
     61def constant_water(t):
     62    return -9.936903
     63
     64#Bo= anuga.shallow_water.boundaries.Transmissive_n_momentum_zero_t_momentum_set_stage_boundary(domain, constant_water) # Outflow for steady uniform flow with a depth integrated velocity of 0.1
     65Bo = anuga.Transmissive_boundary(domain)
    6066domain.set_boundary({'left': Br, 'right': Bo, 'top': Br, 'bottom': Br})
    6167#------------------------------------------------------------------------------
Note: See TracChangeset for help on using the changeset viewer.