Ignore:
Timestamp:
Jun 3, 2009, 1:16:10 PM (16 years ago)
Author:
ole
Message:

Added another memset in quantity_ext.c and fixed some indentation issues to
comply with the Emacs c-style.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • anuga_core/source/anuga/shallow_water/shallow_water_ext.c

    r7143 r7154  
    6868    if (dq1>=dq2){
    6969      if (dq1>=0.0)
    70     *qmax=dq0+dq1;
     70        *qmax=dq0+dq1;
    7171      else
    72     *qmax=dq0;
    73    
     72        *qmax=dq0;
     73     
    7474      *qmin=dq0+dq2;
    7575      if (*qmin>=0.0) *qmin = 0.0;
    … …  
    7777    else{// dq1<dq2
    7878      if (dq2>0)
    79     *qmax=dq0+dq2;
     79        *qmax=dq0+dq2;
    8080      else
    81     *qmax=dq0;
    82    
     81        *qmax=dq0;
     82     
    8383      *qmin=dq0+dq1;   
    8484      if (*qmin>=0.0) *qmin=0.0;
    … …  
    8888    if (dq1<=dq2){
    8989      if (dq1<0.0)
    90     *qmin=dq0+dq1;
     90        *qmin=dq0+dq1;
    9191      else
    92     *qmin=dq0;
    93    
     92        *qmin=dq0;
     93     
    9494      *qmax=dq0+dq2;   
    9595      if (*qmax<=0.0) *qmax=0.0;
    … …  
    9797    else{// dq1>dq2
    9898      if (dq2<0.0)
    99     *qmin=dq0+dq2;
     99        *qmin=dq0+dq2;
    100100      else
    101     *qmin=dq0;
    102    
     101        *qmin=dq0;
     102     
    103103      *qmax=dq0+dq1;
    104104      if (*qmax<=0.0) *qmax=0.0;
    … …  
    683683      // FIXME: Try with this one precomputed
    684684      for (i=0; i<3; i++) {
    685     dz = max(dz, fabs(zv[k3+i]-zc[k]));
     685        dz = max(dz, fabs(zv[k3+i]-zc[k]));
    686686      }
    687687    }
    … …  
    711711     
    712712      if (dz > 0.0) {
    713     alpha = max( min( alpha_balance*hmin/dz, 1.0), 0.0 );     
     713        alpha = max( min( alpha_balance*hmin/dz, 1.0), 0.0 );     
    714714      } else {
    715     alpha = 1.0;  // Flat bed
     715        alpha = 1.0;  // Flat bed
    716716      }
    717717      //printf("Using old style limiter\n");
    … …  
    726726   
    727727      if (hmin < H0) {
    728     alpha = 1.0;
    729     for (i=0; i<3; i++) {
    730      
    731       h_diff = hc_k - hv[i];     
    732       if (h_diff <= 0) {
    733         // Deep water triangle is further away from bed than
    734         // shallow water (hbar < h). Any alpha will do
    735      
    736       } else { 
    737         // Denominator is positive which means that we need some of the
    738         // h-limited stage.
    739        
    740         alpha = min(alpha, (hc_k - H0)/h_diff);     
     728        alpha = 1.0;
     729        for (i=0; i<3; i++) {
     730         
     731          h_diff = hc_k - hv[i];     
     732          if (h_diff <= 0) {
     733            // Deep water triangle is further away from bed than
     734            // shallow water (hbar < h). Any alpha will do
     735     
     736          } else { 
     737            // Denominator is positive which means that we need some of the
     738            // h-limited stage.
     739           
     740            alpha = min(alpha, (hc_k - H0)/h_diff);     
     741          }
     742        }
     743
     744        // Ensure alpha in [0,1]
     745        if (alpha>1.0) alpha=1.0;
     746        if (alpha<0.0) alpha=0.0;
     747       
     748      } else {
     749        // Use w-limited stage exclusively in deeper water.
     750        alpha = 1.0;       
    741751      }
    742752    }
    743 
    744     // Ensure alpha in [0,1]
    745     if (alpha>1.0) alpha=1.0;
    746     if (alpha<0.0) alpha=0.0;
    747    
    748       } else {
    749     // Use w-limited stage exclusively in deeper water.
    750     alpha = 1.0;       
    751       }
    752     }
    753 
    754 
     753   
     754   
    755755    //  Let
    756756    //
    … …  
    770770    //   Momentum is balanced between constant and limited
    771771
    772 
     772   
    773773    if (alpha < 1) {     
    774774      for (i=0; i<3; i++) {
    775      
    776     wv[k3+i] = zv[k3+i] + (1-alpha)*hc_k + alpha*hv[i];
    777 
    778     // Update momentum at vertices
    779     if (use_centroid_velocities == 1) {
    780       // This is a simple, efficient and robust option
    781       // It uses first order approximation of velocities, but retains
    782       // the order used by stage.
    783    
    784       // Speeds at centroids
    785       if (hc_k > epsilon) {
    786         uc = xmomc[k]/hc_k;
    787         vc = ymomc[k]/hc_k;
    788       } else {
    789         uc = 0.0;
    790         vc = 0.0;
    791       }
    792      
    793       // Vertex momenta guaranteed to be consistent with depth guaranteeing
    794       // controlled speed
    795       hv[i] = wv[k3+i] - zv[k3+i]; // Recompute (balanced) vertex depth
    796       xmomv[k3+i] = uc*hv[i];
    797       ymomv[k3+i] = vc*hv[i];
    798      
    799     } else {
    800       // Update momentum as a linear combination of
    801       // xmomc and ymomc (shallow) and momentum
    802       // from extrapolator xmomv and ymomv (deep).
    803       // This assumes that values from xmomv and ymomv have
    804       // been established e.g. by the gradient limiter.
    805 
    806       // FIXME (Ole): I think this should be used with vertex momenta
    807       // computed above using centroid_velocities instead of xmomc
    808       // and ymomc as they'll be more representative first order
    809       // values.
    810      
    811       xmomv[k3+i] = (1-alpha)*xmomc[k] + alpha*xmomv[k3+i];
    812       ymomv[k3+i] = (1-alpha)*ymomc[k] + alpha*ymomv[k3+i];
    813    
    814     }
     775       
     776        wv[k3+i] = zv[k3+i] + (1-alpha)*hc_k + alpha*hv[i];
     777
     778        // Update momentum at vertices
     779        if (use_centroid_velocities == 1) {
     780          // This is a simple, efficient and robust option
     781          // It uses first order approximation of velocities, but retains
     782          // the order used by stage.
     783   
     784          // Speeds at centroids
     785          if (hc_k > epsilon) {
     786            uc = xmomc[k]/hc_k;
     787            vc = ymomc[k]/hc_k;
     788          } else {
     789            uc = 0.0;
     790            vc = 0.0;
     791          }
     792         
     793          // Vertex momenta guaranteed to be consistent with depth guaranteeing
     794          // controlled speed
     795          hv[i] = wv[k3+i] - zv[k3+i]; // Recompute (balanced) vertex depth
     796          xmomv[k3+i] = uc*hv[i];
     797          ymomv[k3+i] = vc*hv[i];
     798     
     799        } else {
     800          // Update momentum as a linear combination of
     801          // xmomc and ymomc (shallow) and momentum
     802          // from extrapolator xmomv and ymomv (deep).
     803          // This assumes that values from xmomv and ymomv have
     804          // been established e.g. by the gradient limiter.
     805         
     806          // FIXME (Ole): I think this should be used with vertex momenta
     807          // computed above using centroid_velocities instead of xmomc
     808          // and ymomc as they'll be more representative first order
     809          // values.
     810         
     811          xmomv[k3+i] = (1-alpha)*xmomc[k] + alpha*xmomv[k3+i];
     812          ymomv[k3+i] = (1-alpha)*ymomc[k] + alpha*ymomv[k3+i];
     813         
     814        }
    815815      }
    816816    }
    … …  
    843843      if (hc < minimum_allowed_height) {
    844844       
    845     // Set momentum to zero and ensure h is non negative
    846     xmomc[k] = 0.0;
    847     ymomc[k] = 0.0;
    848     if (hc <= 0.0) wc[k] = zc[k];
     845        // Set momentum to zero and ensure h is non negative
     846        xmomc[k] = 0.0;
     847        ymomc[k] = 0.0;
     848        if (hc <= 0.0) wc[k] = zc[k];
    849849      }
    850850    }
    … …  
    854854    for (k=0; k<N; k++) {
    855855      hc = wc[k] - zc[k];
    856 
     856     
    857857      if (hc < minimum_allowed_height) {
    858858
    … …  
    866866        } else {
    867867          //Reduce excessive speeds derived from division by small hc
    868         //FIXME (Ole): This may be unnecessary with new slope limiters
    869         //in effect.
     868          //FIXME (Ole): This may be unnecessary with new slope limiters
     869          //in effect.
    870870         
    871871          u = xmomc[k]/hc;
    872       if (fabs(u) > maximum_allowed_speed) {
    873         reduced_speed = maximum_allowed_speed * u/fabs(u);
    874         //printf("Speed (u) has been reduced from %.3f to %.3f\n",
    875         //   u, reduced_speed);
    876         xmomc[k] = reduced_speed * hc;
    877       }
    878 
     872          if (fabs(u) > maximum_allowed_speed) {
     873            reduced_speed = maximum_allowed_speed * u/fabs(u);
     874            //printf("Speed (u) has been reduced from %.3f to %.3f\n",
     875            //   u, reduced_speed);
     876            xmomc[k] = reduced_speed * hc;
     877          }
     878         
    879879          v = ymomc[k]/hc;
    880       if (fabs(v) > maximum_allowed_speed) {
    881         reduced_speed = maximum_allowed_speed * v/fabs(v);
    882         //printf("Speed (v) has been reduced from %.3f to %.3f\n",
    883         //   v, reduced_speed);
    884         ymomc[k] = reduced_speed * hc;
    885       }
     880          if (fabs(v) > maximum_allowed_speed) {
     881            reduced_speed = maximum_allowed_speed * v/fabs(v);
     882            //printf("Speed (v) has been reduced from %.3f to %.3f\n",
     883            //   v, reduced_speed);
     884            ymomc[k] = reduced_speed * hc;
     885          }
    886886        }
    887887      }
Note: See TracChangeset for help on using the changeset viewer.