Changeset 4251


Ignore:
Timestamp:
Feb 9, 2007, 4:23:15 PM (18 years ago)
Author:
ole
Message:

Pulled speed computation out into a separate function that can be used by all flux functions.

File:
1 edited

Legend:

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

    r4239 r4251  
    125125}
    126126
     127
     128// Function to obtain speed from momentum and depth.
     129// This is used by flux functions
     130// Input parameters uh and h may be modified by this function.
     131double _compute_speed(double *uh,
     132                      double *h,
     133                      double epsilon,
     134                      double h0) {
     135 
     136  double u;
     137 
     138  if (*h < epsilon) {
     139    *h = 0.0;  //Could have been negative
     140    u = 0.0;
     141  } else {
     142    u = *uh/(*h + h0/ *h);
     143  }
     144 
     145  return u;
     146}
     147
    127148// Computational function for flux computation (using stage w=z+h)
    128149int flux_function_central(double *q_left, double *q_right,
     
    167188
    168189  //Compute speeds in x-direction
    169   w_left = q_left_copy[0];              // h+z
     190  w_left = q_left_copy[0];         
    170191  h_left = w_left-z;
    171192  uh_left = q_left_copy[1];
    172 
    173   if (h_left < epsilon) {
    174     h_left = 0.0;  //Could have been negative
    175     u_left = 0.0;
    176   } else {
    177     u_left = uh_left/(h_left + h0/h_left);
    178   }
     193  u_left =_compute_speed(&uh_left, &h_left, epsilon, h0);
    179194
    180195  w_right = q_right_copy[0];
    181196  h_right = w_right-z;
    182197  uh_right = q_right_copy[1];
    183 
    184   if (h_right < epsilon) {
    185     h_right = 0.0; //Could have been negative
    186     u_right = 0.0;
    187   } else {
    188     u_right = uh_right/(h_right + h0/h_right);
    189   }
     198  u_right =_compute_speed(&uh_right, &h_right, epsilon, h0); 
    190199
    191200  //Momentum in y-direction
     
    288297
    289298  //Compute speeds in x-direction
    290   w_left = q_left_copy[0];              // h+z
     299  w_left = q_left_copy[0];         
    291300  h_left = w_left-z;
    292301  uh_left = q_left_copy[1];
    293 
    294   if (h_left < epsilon) {
    295     h_left = 0.0;  //Could have been negative
    296     u_left = 0.0;
    297   } else {
    298     u_left = uh_left/(h_left + h0/h_left);
    299   }
     302  u_left =_compute_speed(&uh_left, &h_left, epsilon, h0);
    300303
    301304  w_right = q_right_copy[0];
    302305  h_right = w_right-z;
    303306  uh_right = q_right_copy[1];
    304 
    305   if (h_right < epsilon) {
    306     h_right = 0.0; //Could have been negative
    307     u_right = 0.0;
    308   } else {
    309     u_right = uh_right/(h_right + h0/h_right);
    310   }
     307  u_right =_compute_speed(&uh_right, &h_right, epsilon, h0); 
     308
    311309
    312310  //Momentum in y-direction
Note: See TracChangeset for help on using the changeset viewer.