Changeset 5297 for anuga_core/source


Ignore:
Timestamp:
May 9, 2008, 9:44:37 AM (16 years ago)
Author:
ole
Message:

Implemented option use_centroid_velocities as per ticket:262.
The run_profile script on a linux amd 64 (bogong) improved from 37s to
33s or about 10%.

File:
1 edited

Legend:

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

    r5290 r5297  
    550550                              double H0,
    551551                              int tight_slope_limiters,
    552                               int use_centroid_velocities,                           
     552                              int use_centroid_velocities,                   
    553553                              double alpha_balance) {
    554554
     
    652652      }
    653653
    654       //printf("alpha=%f, tight_slope_limiters=%d\n", alpha, tight_slope_limiters);
    655 
    656       /*     
    657       // Experimental code for controlling velocities at vertices.
    658       // Adjust alpha (down towards first order) such that
    659       // velocities at vertices remain within one order of magnitude
    660       // of those at the centroid.
    661      
    662       for (i=0; i<3; i++) {     
    663      
    664         // FIXME (Ole): Simplify (remove) when (if) hvbar gets retired
    665         if (beta_h > epsilon) {
    666           hc_k = hvbar[k3+i]; // Depth to be used at vertices
    667         }
    668        
    669        
    670         alpha = velocity_balance(xmomv[k3+i], xmomc[k],
    671                                  hv[i], hc_k, alpha, epsilon);
    672                                  
    673         alpha = velocity_balance(ymomv[k3+i], ymomc[k],
    674                                  hv[i], hc_k, alpha, epsilon);                           
    675       }
    676       */
    677654    }
    678655           
    679     //printf("k=%d, hmin=%.2f, dz=%.2f, alpha=%.2f, alpha_balance=%.2f\n",
    680     //     k, hmin, dz, alpha, alpha_balance);
    681 
    682     //printf("dz = %.3f, alpha = %.8f\n", dz, alpha);
    683656
    684657    //  Let
     
    711684
    712685        // Update momentum at vertices
    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.
    717        
    718           // Speeds at centroids
    719           if (hc_k > epsilon) {
    720             uc = xmomc[k]/hc_k;
    721             vc = ymomc[k]/hc_k;
    722           } else {
    723             uc = 0.0;
    724             vc = 0.0;
    725           }
    726          
    727           // Vertex momenta guaranteed to be consistent with depth guaranteeing
    728           // controlled speed
    729           hv[i] = wv[k3+i] - zv[k3+i]; // Recompute (balanced) vertex depth
    730           xmomv[k3+i] = uc*hv[i];
    731           ymomv[k3+i] = vc*hv[i];
    732          
    733         } else {
     686        if (use_centroid_velocities == 0) {
    734687          // Update momentum as a linear combination of
    735688          // xmomc and ymomc (shallow) and momentum
     
    744697        }
    745698      }
     699    } // If alpha == 1 use quantities as calculated by the gradient-limiter
     700   
     701    if (use_centroid_velocities == 1) {   
     702      // This is a simple, efficient and robust option in shallow water
     703      // It uses first order approximation of velocities, but retains
     704      // the order used by stage.
     705      // In this case the xmomv and ymomv calculations should be switched off
     706      // in the gradient limiter.
     707
     708      for (i=0; i<3; i++) {   
     709       
     710        // Speeds at centroids
     711        if (hc_k > epsilon) {
     712          uc = xmomc[k]/hc_k;
     713          vc = ymomc[k]/hc_k;
     714        } else {
     715          uc = 0.0;
     716          vc = 0.0;
     717        }
     718       
     719        // Vertex momenta guaranteed to be consistent with depth guaranteeing
     720        // controlled speed
     721        hv[i] = wv[k3+i] - zv[k3+i]; // Recompute (balanced) vertex depth
     722        xmomv[k3+i] = uc*hv[i];
     723        ymomv[k3+i] = vc*hv[i];
     724      }
    746725    }
    747726  }
     
    10481027  // Local variables
    10491028  double a, b; // Gradient vector used to calculate vertex values from centroids
    1050   int k,k0,k1,k2,k3,k6,coord_index,i;
     1029  int k,k0,k1,k2,k3,k6,coord_index, i;
    10511030  double x,y,x0,y0,x1,y1,x2,y2,xv0,yv0,xv1,yv1,xv2,yv2; // Vertices of the auxiliary triangle
    10521031  double dx1,dx2,dy1,dy2,dxv0,dxv1,dxv2,dyv0,dyv1,dyv2,dq0,dq1,dq2,area2;
     
    12131192      for (i=0;i<3;i++)
    12141193        stage_vertex_values[k3+i]=stage_centroid_values[k]+dqv[i];
     1194     
     1195     
     1196      if (use_centroid_velocities == 1) {
     1197        // Use first order reconstruction using speeds only.
     1198       
     1199        // This happens in balance_deep_and_shallow so there
     1200        // is no need to do more here
     1201       
     1202        // FIXME (Ole): Optionally we could put the computation here but then
     1203        // it'd go into all other variants of the gradient limiter 
     1204        continue;
     1205      }
    12151206     
    12161207     
Note: See TracChangeset for help on using the changeset viewer.