Ignore:
Timestamp:
May 7, 2008, 5:40:11 PM (16 years ago)
Author:
ole
Message:

Refactored options, so that the use of centroid_velocities is a separate
option to tight_slope_limiters.

File:
1 edited

Legend:

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

    r5238 r5290  
    144144                         
    145145  // Compute Froude number; v/sqrt(gh)
    146  
     146  // FIXME (Ole): Not currently in use
     147   
    147148  double froude_number;
    148149 
     
    549550                              double H0,
    550551                              int tight_slope_limiters,
     552                              int use_centroid_velocities,                           
    551553                              double alpha_balance) {
    552554
    553   int k, k3, i, excessive_froude_number=0;
     555  int k, k3, i; //, excessive_froude_number=0;
    554556
    555557  double dz, hmin, alpha, h_diff, hc_k;
    556558  double epsilon = 1.0e-6; // FIXME: Temporary measure
    557   double g = 9.81; // FIXME: Temporary measure
    558   double hv[3], h; // Depths at vertices
    559   double Fx, Fy; // Froude numbers
     559  double hv[3]; // Depths at vertices
     560  //double Fx, Fy; // Froude numbers
    560561  double uc, vc; // Centroid speeds
    561562
     
    701702    if (alpha < 1) {     
    702703      for (i=0; i<3; i++) {
     704     
    703705        // FIXME (Ole): Simplify when (if) hvbar gets retired       
    704706        if (beta_h > epsilon) {   
     
    709711
    710712        // Update momentum at vertices
    711 
    712 
    713         // Update momentum at vertices
    714         // Update momentum as a linear combination of
    715         // xmomc and ymomc (shallow) and momentum
    716         // from extrapolator xmomv and ymomv (deep).
    717 
    718 
    719         //xmomv[k3+i] = (1-alpha)*xmomc[k] + alpha*xmomv[k3+i];
    720         //ymomv[k3+i] = (1-alpha)*ymomc[k] + alpha*ymomv[k3+i];
    721 
    722         if (tight_slope_limiters == 1) {
    723           // FIXME(Ole): Here's what I think (as of 17 Nov 2007)
    724           // we need to do. Simple and efficient:
     713        if (use_centroid_velocities == 1) {
     714          // This is a simple, efficient and robust option
     715          // It uses first order approximation of velocities, but retains
     716          // the order used by stage.
    725717       
    726718          // Speeds at centroids
     
    738730          xmomv[k3+i] = uc*hv[i];
    739731          ymomv[k3+i] = vc*hv[i];
     732         
    740733        } else {
    741734          // Update momentum as a linear combination of
    742735          // xmomc and ymomc (shallow) and momentum
    743736          // from extrapolator xmomv and ymomv (deep).
    744           // FIXME (Ole): Is this really needed? Could we use the above
    745           // instead?
     737          // This assumes that values from xmomv and ymomv have
     738          // been established e.g. by the gradient limiter.
     739
    746740         
    747741          xmomv[k3+i] = (1-alpha)*xmomc[k] + alpha*xmomv[k3+i];
     
    749743       
    750744        }
    751       }
    752     }
    753        
    754 
    755        
    756     if (0) {  // FIXME(Ole): Disabled while testing balancing of velocities above
    757       //if (tight_slope_limiters == 1) {               
    758    
    759       // Ensure that the Froude number is kept realistic at vertices
    760       // FIXME (Ole): I think it could be used to adjust alpha down
    761       // whenever Fr is too large. Possible make sure it doesn't deviate
    762       // too much from the value at the centroid. I like this idea!
    763      
    764       // FIXME (Ole): currently only used with tights_SL
    765      
    766       // FIXME (Ole): may not be necessary now
    767 
    768       excessive_froude_number=0;   
    769       for (i=0; i<3; i++) {   
    770      
    771         // Recalculate depth at vertex i
    772         h = wv[k3+i] - zv[k3+i];
    773        
    774         Fx = compute_froude_number(xmomv[k3+i], h, g, epsilon);
    775         Fy = compute_froude_number(ymomv[k3+i], h, g, epsilon);
    776        
    777         if ( (fabs(Fx) > 100.0) || (fabs(Fy) > 100.0)) {
    778           // FIXME: use max_froude - or base it on centroid value of F
    779           // printf("Excessive Froude number detected: %f or %f\n", Fx, Fy);
    780           excessive_froude_number=1;         
    781         }
    782       }
    783      
    784       if (excessive_froude_number) {
    785    
    786         // printf("Adjusting momentum to first order.\n");
    787         // Go back to first order (alpha = 0) for this triangle
    788         for (i=0; i<3; i++) {         
    789           xmomv[k3+i] = xmomc[k];
    790           ymomv[k3+i] = ymomc[k];
    791          
    792           if (beta_h > epsilon) {         
    793             wv[k3+i] = zv[k3+i] + hvbar[k3+i];
    794           } else {
    795             wv[k3+i] = zv[k3+i] + hc_k;
    796           }
    797         }       
    798745      }
    799746    }
     
    10941041                                  double* ymom_vertex_values,
    10951042                                  double* elevation_vertex_values,
    1096                                   int optimise_dry_cells) {
     1043                                  int optimise_dry_cells,
     1044                                  int use_centroid_velocities) {
    10971045                                 
    10981046                                 
     
    13591307      // One internal neighbour and gradient is in direction of the neighbour's centroid
    13601308     
    1361       // Find the only internal neighbour
     1309      // Find the only internal neighbour (k1?)
    13621310      for (k2=k3;k2<k3+3;k2++){
    13631311        // Find internal neighbour of triangle k     
     
    13671315          break;
    13681316      }
     1317     
    13691318      if ((k2==k3+3)) {
    13701319        // If we didn't find an internal neighbour
     
    13961345      dy2=dx2*dy1;
    13971346      dx2*=dx1;
    1398      
    13991347     
    14001348     
     
    15611509  double beta_w, beta_w_dry, beta_uh, beta_uh_dry, beta_vh, beta_vh_dry;   
    15621510  double minimum_allowed_height, epsilon;
    1563   int optimise_dry_cells, number_of_elements, e;
     1511  int optimise_dry_cells, number_of_elements, e, use_centroid_velocities;
    15641512 
    15651513  // Provisional jumps from centroids to v'tices and safety factor re limiting
    15661514  // by which these jumps are limited
    15671515  // Convert Python arguments to C
    1568   if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOi",
     1516  if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOii",
    15691517                        &domain,
    15701518                        &surrogate_neighbours,
     
    15801528                        &ymom_vertex_values,
    15811529                        &elevation_vertex_values,
    1582                         &optimise_dry_cells)) {                 
     1530                        &optimise_dry_cells,
     1531                        &use_centroid_velocities)) {                   
    15831532                       
    15841533    PyErr_SetString(PyExc_RuntimeError,
     
    16251574                                   (double*) ymom_vertex_values -> data,
    16261575                                   (double*) elevation_vertex_values -> data,
    1627                                    optimise_dry_cells);
     1576                                   optimise_dry_cells,
     1577                                   use_centroid_velocities);
    16281578  if (e == -1) {
    16291579    // Use error string set inside computational routine
     
    22712221  double H0, beta_h;
    22722222
    2273   int N, tight_slope_limiters; //, err;
     2223  int N, tight_slope_limiters, use_centroid_velocities; //, err;
    22742224
    22752225  // Convert Python arguments to C
     
    23182268 
    23192269 
     2270  Tmp = PyObject_GetAttrString(domain, "use_centroid_velocities");
     2271  if (!Tmp) {
     2272    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not obtain object use_centroid_velocities from domain");
     2273    return NULL;
     2274  } 
     2275  use_centroid_velocities = PyInt_AsLong(Tmp);
     2276  Py_DECREF(Tmp);
     2277 
     2278 
     2279 
    23202280  N = wc -> dimensions[0];
    23212281  _balance_deep_and_shallow(N,
     
    23322292                            H0,
    23332293                            (int) tight_slope_limiters,
     2294                            (int) use_centroid_velocities,                         
    23342295                            alpha_balance);
    23352296
Note: See TracChangeset for help on using the changeset viewer.