Changeset 4823


Ignore:
Timestamp:
Nov 16, 2007, 5:17:49 PM (16 years ago)
Author:
ole
Message:

Work towards tight_slope_limiters with balanced velocities and/or Froude numbers.

Location:
anuga_core/source/anuga/shallow_water
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/shallow_water/shallow_water_ext.c

    r4815 r4823  
    464464*/
    465465
     466
     467
     468
     469double velocity_balance(double uh_i,
     470                        double uh,
     471                        double h_i,
     472                        double h,
     473                        double alpha,
     474                        double epsilon) {
     475  // Find alpha such that speed at the vertex is within one
     476  // order of magnitude of the centroid speed   
     477
     478  // FIXME(Ole): Work in progress
     479 
     480  double a, b, estimate;
     481  double m=10; // One order of magnitude - allow for velocity deviations at vertices
     482
     483 
     484  printf("alpha = %f, uh_i=%f, uh=%f, h_i=%f, h=%f\n",
     485         alpha, uh_i, uh, h_i, h);
     486     
     487 
     488 
     489   
     490  // Shorthands and determine inequality
     491  if (fabs(uh) < epsilon) {
     492    a = 1.0e10; // Limit
     493  } else {
     494    a = fabs(uh_i - uh)/fabs(uh);
     495  }
     496     
     497  if (h < epsilon) {
     498    b = 1.0e10; // Limit
     499  } else {       
     500    b = m*fabs(h_i - h)/h;
     501  }
     502
     503  printf("a %f, b %f\n", a, b);
     504
     505  if (a > b) {
     506    estimate = (m-1)/(a-b);             
     507   
     508    printf("Alpha %f, estimate %f\n",
     509           alpha, estimate);   
     510           
     511    if (alpha < estimate) {
     512      printf("Adjusting alpha from %f to %f\n",
     513             alpha, estimate);
     514      alpha = estimate;
     515    }
     516  } else {
     517 
     518    if (h < h_i) {
     519      estimate = (m-1)/(a-b);                 
     520   
     521      printf("Alpha %f, estimate %f\n",
     522             alpha, estimate);   
     523           
     524      if (alpha < estimate) {
     525        printf("Adjusting alpha from %f to %f\n",
     526               alpha, estimate);
     527        alpha = estimate;
     528      }   
     529    }
     530    // Always fulfilled as alpha and m-1 are non negative
     531  }
     532 
     533 
     534  return alpha;
     535}
     536
     537
    466538int _balance_deep_and_shallow(int N,
    467539                              double beta_h,
     
    480552
    481553  int k, k3, i, excessive_froude_number=0;
     554
    482555  double dz, hmin, alpha, h_diff, hc_k;
    483556  double epsilon = 1.0e-6; // FIXME: Temporary measure
     
    539612   
    540613      // Make alpha as large as possible but still ensure that
    541       // final depth is positive
     614      // final depth is positive and that velocities at vertices
     615      // are controlled
    542616   
    543617      if (hmin < H0) {
    544618        alpha = 1.0;
    545619        for (i=0; i<3; i++) {
    546 
    547           // FIXME (Ole): Simplify when (if) hvbar gets retired
     620       
     621          // FIXME (Ole): Simplify (remove) when (if) hvbar gets retired
    548622          if (beta_h > epsilon) {
    549             h_diff = hvbar[k3+i] - hv[i];
    550           } else {
    551             h_diff = hc_k - hv[i];       
     623            hc_k = hvbar[k3+i]; // Depth to be used at vertices
    552624          }
    553        
     625         
     626          h_diff = hc_k - hv[i];         
    554627          if (h_diff <= 0) {
    555628            // Deep water triangle is further away from bed than
     
    560633            // h-limited stage.
    561634           
    562             // FIXME (Ole): Simplify when (if) hvbar gets retired           
    563             if (beta_h > epsilon) {       
    564               alpha = min(alpha, (hvbar[k3+i] - H0)/h_diff);
    565             } else {
    566               alpha = min(alpha, (hc_k - H0)/h_diff);       
    567             }
     635            alpha = min(alpha, (hc_k - H0)/h_diff);         
    568636          }
    569637        }
     
    574642       
    575643      } else {
    576         // Use w-limited stage exclusively
     644        // Use w-limited stage exclusively in deeper water.
    577645        alpha = 1.0;       
    578646      }
     647
     648     
     649      /*     
     650      // Experimental code for controlling velocities at vertices.
     651      // Adjust alpha (down towards first order) such that
     652      // velocities at vertices remain within one order of magnitude
     653      // of those at the centroid.
     654     
     655      for (i=0; i<3; i++) {     
     656     
     657        // FIXME (Ole): Simplify (remove) when (if) hvbar gets retired
     658        if (beta_h > epsilon) {
     659          hc_k = hvbar[k3+i]; // Depth to be used at vertices
     660        }
     661       
     662       
     663        alpha = velocity_balance(xmomv[k3+i], xmomc[k],
     664                                 hv[i], hc_k, alpha, epsilon);
     665                                 
     666        alpha = velocity_balance(ymomv[k3+i], ymomc[k],
     667                                 hv[i], hc_k, alpha, epsilon);                           
     668      }
     669      */
    579670    }
    580        
    581    
    582        
     671           
    583672    //printf("k=%d, hmin=%.2f, dz=%.2f, alpha=%.2f, alpha_balance=%.2f\n",
    584673    //     k, hmin, dz, alpha, alpha_balance);
     
    624713
    625714       
     715    //if (0) {  // FIXME(Ole): Disabled while testing balancing of velocities above
    626716    if (tight_slope_limiters == 1) {                   
    627717   
  • anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

    r4815 r4823  
    28612861        u = xmomentum/depth
    28622862        v = ymomentum/depth
     2863
     2864        # Assert that all vertex velocities stay within one
     2865        # order of magnitude of centroid velocities.
     2866        #print u.vertex_values[1,:]
     2867        #print u.centroid_values[1]
     2868       
     2869        assert alltrue(absolute(u.vertex_values[1,:]) <= absolute(u.centroid_values[1])*10)
     2870        assert alltrue(absolute(v.vertex_values[1,:]) <= absolute(v.centroid_values[1])*10)
    28632871
    28642872        denom = (depth*g)**0.5
Note: See TracChangeset for help on using the changeset viewer.