Ignore:
Timestamp:
Apr 24, 2009, 5:22:14 PM (16 years ago)
Author:
rwilson
Message:

Back-merge from Numeric trunk to numpy branch.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/numpy/anuga/shallow_water/shallow_water_ext.c

    r6780 r6902  
    2828
    2929// Computational function for rotation
     30// FIXME: Perhaps inline this and profile
    3031int _rotate(double *q, double n1, double n2) {
    3132  /*Rotate the momentum component q (q[1], q[2])
    … …  
    5253
    5354int find_qmin_and_qmax(double dq0, double dq1, double dq2,
    54                        double *qmin, double *qmax){
     55               double *qmin, double *qmax){
    5556  // Considering the centroid of an FV triangle and the vertices of its
    5657  // auxiliary triangle, find
    … …  
    6768    if (dq1>=dq2){
    6869      if (dq1>=0.0)
    69         *qmax=dq0+dq1;
     70    *qmax=dq0+dq1;
    7071      else
    71         *qmax=dq0;
    72        
     72    *qmax=dq0;
     73   
    7374      *qmin=dq0+dq2;
    7475      if (*qmin>=0.0) *qmin = 0.0;
    … …  
    7677    else{// dq1<dq2
    7778      if (dq2>0)
    78         *qmax=dq0+dq2;
     79    *qmax=dq0+dq2;
    7980      else
    80         *qmax=dq0;
    81        
    82       *qmin=dq0+dq1;   
     81    *qmax=dq0;
     82   
     83      *qmin=dq0+dq1;   
    8384      if (*qmin>=0.0) *qmin=0.0;
    8485    }
    … …  
    8788    if (dq1<=dq2){
    8889      if (dq1<0.0)
    89         *qmin=dq0+dq1;
     90    *qmin=dq0+dq1;
    9091      else
    91         *qmin=dq0;
    92        
    93       *qmax=dq0+dq2;   
     92    *qmin=dq0;
     93   
     94      *qmax=dq0+dq2;   
    9495      if (*qmax<=0.0) *qmax=0.0;
    9596    }
    9697    else{// dq1>dq2
    9798      if (dq2<0.0)
    98         *qmin=dq0+dq2;
     99    *qmin=dq0+dq2;
    99100      else
    100         *qmin=dq0;
    101        
     101    *qmin=dq0;
     102   
    102103      *qmax=dq0+dq1;
    103104      if (*qmax<=0.0) *qmax=0.0;
    … …  
    133134 
    134135  phi=min(r*beta_w,1.0);
    135   for (i=0;i<3;i++)
    136     dqv[i]=dqv[i]*phi;
     136  //for (i=0;i<3;i++)
     137  dqv[0]=dqv[0]*phi;
     138  dqv[1]=dqv[1]*phi;
     139  dqv[2]=dqv[2]*phi;
    137140   
    138141  return 0;
    … …  
    141144
    142145double compute_froude_number(double uh,
    143                              double h,
    144                              double g,
    145                              double epsilon) {
    146                          
     146                 double h,
     147                 double g,
     148                 double epsilon) {
     149             
    147150  // Compute Froude number; v/sqrt(gh)
    148151  // FIXME (Ole): Not currently in use
    … …  
    166169// This is used by flux functions
    167170// Input parameters uh and h may be modified by this function.
     171
     172// FIXME: Perhaps inline this and profile
    168173double _compute_speed(double *uh,
    169                       double *h,
    170                       double epsilon,
    171                       double h0) {
     174              double *h,
     175              double epsilon,
     176              double h0) {
    172177 
    173178  double u;
    … …  
    188193}
    189194
     195
     196// Optimised squareroot computation
     197//
     198//See
     199//http://www.lomont.org/Math/Papers/2003/InvSqrt.pdf
     200//and http://mail.opensolaris.org/pipermail/tools-gcc/2005-August/000047.html
     201float fast_squareroot_approximation(float number) {
     202  float x;
     203  const float f = 1.5;
     204
     205  // Allow i and y to occupy the same memory
     206  union
     207  {
     208    long i;
     209    float y;
     210  } u;
     211 
     212  // Good initial guess 
     213  x = number * 0.5; 
     214  u.y  = number;
     215  u.i  = 0x5f3759df - ( u.i >> 1 );
     216 
     217  // Take a few iterations
     218  u.y  = u.y * ( f - ( x * u.y * u.y ) );
     219  u.y  = u.y * ( f - ( x * u.y * u.y ) );
     220   
     221  return number * u.y;
     222}
     223
     224
     225
     226// Optimised squareroot computation (double version, slower)
     227double Xfast_squareroot_approximation(double number) {
     228  double x;
     229  const double f = 1.5;
     230   
     231  // Allow i and y to occupy the same memory   
     232  union
     233  {
     234    long long i;
     235    double y;
     236  } u;
     237
     238  // Good initial guess
     239  x = number * 0.5;
     240  u.y  = number;
     241  u.i  = 0x5fe6ec85e7de30daLL - ( u.i >> 1 );
     242 
     243  // Take a few iterations 
     244  u.y  = u.y * ( f - ( x * u.y * u.y ) );
     245  u.y  = u.y * ( f - ( x * u.y * u.y ) );
     246
     247  return number * u.y;
     248}
     249
     250
    190251// Innermost flux function (using stage w=z+h)
    191252int _flux_function_central(double *q_left, double *q_right,
    192                            double z_left, double z_right,
    193                            double n1, double n2,
    194                            double epsilon, double H0, double g,
    195                            double *edgeflux, double *max_speed) {
     253                           double z_left, double z_right,
     254                           double n1, double n2,
     255                           double epsilon, double H0, double g,
     256                           double *edgeflux, double *max_speed)
     257{
    196258
    197259  /*Compute fluxes between volumes for the shallow water wave equation
    … …  
    212274  double v_left, v_right; 
    213275  double s_min, s_max, soundspeed_left, soundspeed_right;
    214   double denom, z;
     276  double denom, inverse_denominator, z;
    215277
    216278  // Workspace (allocate once, use many)
    … …  
    219281  double h0 = H0*H0; // This ensures a good balance when h approaches H0.
    220282                     // But evidence suggests that h0 can be as little as
    221                      // epsilon!
     283             // epsilon!
    222284 
    223285  // Copy conserved quantities to protect from modification
    224   for (i=0; i<3; i++) {
    225     q_left_rotated[i] = q_left[i];
    226     q_right_rotated[i] = q_right[i];
    227   }
     286  q_left_rotated[0] = q_left[0];
     287  q_right_rotated[0] = q_right[0];
     288  q_left_rotated[1] = q_left[1];
     289  q_right_rotated[1] = q_right[1];
     290  q_left_rotated[2] = q_left[2];
     291  q_right_rotated[2] = q_right[2];
    228292
    229293  // Align x- and y-momentum with x-axis
    … …  
    231295  _rotate(q_right_rotated, n1, n2);
    232296
    233   z = (z_left+z_right)/2; // Average elevation values.
    234                           // Even though this will nominally allow for discontinuities
    235                           // in the elevation data, there is currently no numerical
    236                           // support for this so results may be strange near jumps in the bed.
     297  z = 0.5*(z_left + z_right); // Average elevation values.
     298                            // Even though this will nominally allow for discontinuities
     299                            // in the elevation data, there is currently no numerical
     300                            // support for this so results may be strange near jumps in the bed.
    237301
    238302  // Compute speeds in x-direction
    239303  w_left = q_left_rotated[0];         
    240   h_left = w_left-z;
     304  h_left = w_left - z;
    241305  uh_left = q_left_rotated[1];
    242306  u_left = _compute_speed(&uh_left, &h_left, epsilon, h0);
    243307
    244308  w_right = q_right_rotated[0];
    245   h_right = w_right-z;
     309  h_right = w_right - z;
    246310  uh_right = q_right_rotated[1];
    247311  u_right = _compute_speed(&uh_right, &h_right, epsilon, h0); 
    … …  
    257321  // Maximal and minimal wave speeds
    258322  soundspeed_left  = sqrt(g*h_left);
    259   soundspeed_right = sqrt(g*h_right);
    260 
    261   s_max = max(u_left+soundspeed_left, u_right+soundspeed_right);
    262   if (s_max < 0.0) s_max = 0.0;
    263 
    264   s_min = min(u_left-soundspeed_left, u_right-soundspeed_right);
    265   if (s_min > 0.0) s_min = 0.0;
    266 
     323  soundspeed_right = sqrt(g*h_right); 
     324 
     325  // Code to use fast square root optimisation if desired.
     326  //soundspeed_left  = fast_squareroot_approximation(g*h_left);
     327  //soundspeed_right = fast_squareroot_approximation(g*h_right);
     328
     329  s_max = max(u_left + soundspeed_left, u_right + soundspeed_right);
     330  if (s_max < 0.0)
     331  {
     332    s_max = 0.0;
     333  }
     334
     335  s_min = min(u_left - soundspeed_left, u_right - soundspeed_right);
     336  if (s_min > 0.0)
     337  {
     338    s_min = 0.0;
     339  }
    267340 
    268341  // Flux formulas
    … …  
    275348  flux_right[2] = u_right*vh_right;
    276349
    277  
    278350  // Flux computation
    279   denom = s_max-s_min;
    280   if (denom < epsilon) { // FIXME (Ole): Try using H0 here
    281     for (i=0; i<3; i++) edgeflux[i] = 0.0;
     351  denom = s_max - s_min;
     352  if (denom < epsilon)
     353  { // FIXME (Ole): Try using H0 here
     354    memset(edgeflux, 0, 3*sizeof(double));
    282355    *max_speed = 0.0;
    283   } else {
    284     for (i=0; i<3; i++) {
     356  }
     357  else
     358  {
     359    inverse_denominator = 1.0/denom;
     360    for (i = 0; i < 3; i++)
     361    {
    285362      edgeflux[i] = s_max*flux_left[i] - s_min*flux_right[i];
    286       edgeflux[i] += s_max*s_min*(q_right_rotated[i]-q_left_rotated[i]);
    287       edgeflux[i] /= denom;
     363      edgeflux[i] += s_max*s_min*(q_right_rotated[i] - q_left_rotated[i]);
     364      edgeflux[i] *= inverse_denominator;
    288365    }
    289366
    … …  
    315392// Computational function for flux computation (using stage w=z+h)
    316393int flux_function_kinetic(double *q_left, double *q_right,
    317                   double z_left, double z_right,
    318                   double n1, double n2,
    319                   double epsilon, double H0, double g,
    320                   double *edgeflux, double *max_speed) {
     394          double z_left, double z_right,
     395          double n1, double n2,
     396          double epsilon, double H0, double g,
     397          double *edgeflux, double *max_speed) {
    321398
    322399  /*Compute fluxes between volumes for the shallow water wave equation
    … …  
    413490
    414491void _manning_friction(double g, double eps, int N,
    415                        double* w, double* z,
    416                        double* uh, double* vh,
    417                        double* eta, double* xmom, double* ymom) {
     492               double* w, double* z,
     493               double* uh, double* vh,
     494               double* eta, double* xmom, double* ymom) {
    418495
    419496  int k;
    … …  
    426503        S = -g * eta[k]*eta[k] * sqrt((uh[k]*uh[k] + vh[k]*vh[k]));
    427504        S /= pow(h, 7.0/3);      //Expensive (on Ole's home computer)
    428         //S /= exp(7.0/3.0*log(h));      //seems to save about 15% over manning_friction
     505        //S /= exp((7.0/3.0)*log(h));      //seems to save about 15% over manning_friction
    429506        //S /= h*h*(1 + h/3.0 - h*h/9.0); //FIXME: Could use a Taylor expansion
    430507
    … …  
    441518/*
    442519void _manning_friction_explicit(double g, double eps, int N,
    443                        double* w, double* z,
    444                        double* uh, double* vh,
    445                        double* eta, double* xmom, double* ymom) {
     520               double* w, double* z,
     521               double* uh, double* vh,
     522               double* eta, double* xmom, double* ymom) {
    446523
    447524  int k;
    … …  
    452529      h = w[k]-z[k];
    453530      if (h >= eps) {
    454         S = -g * eta[k]*eta[k] * sqrt((uh[k]*uh[k] + vh[k]*vh[k]));
    455         S /= pow(h, 7.0/3);      //Expensive (on Ole's home computer)
    456         //S /= exp(7.0/3.0*log(h));      //seems to save about 15% over manning_friction
    457         //S /= h*h*(1 + h/3.0 - h*h/9.0); //FIXME: Could use a Taylor expansion
    458 
    459 
    460         //Update momentum
    461         xmom[k] += S*uh[k];
    462         ymom[k] += S*vh[k];
     531    S = -g * eta[k]*eta[k] * sqrt((uh[k]*uh[k] + vh[k]*vh[k]));
     532    S /= pow(h, 7.0/3);      //Expensive (on Ole's home computer)
     533    //S /= exp(7.0/3.0*log(h));      //seems to save about 15% over manning_friction
     534    //S /= h*h*(1 + h/3.0 - h*h/9.0); //FIXME: Could use a Taylor expansion
     535
     536
     537    //Update momentum
     538    xmom[k] += S*uh[k];
     539    ymom[k] += S*vh[k];
    463540      }
    464541    }
    … …  
    471548
    472549double velocity_balance(double uh_i,
    473                         double uh,
    474                         double h_i,
    475                         double h,
    476                         double alpha,
    477                         double epsilon) {
     550            double uh,
     551            double h_i,
     552            double h,
     553            double alpha,
     554            double epsilon) {
    478555  // Find alpha such that speed at the vertex is within one
    479   // order of magnitude of the centroid speed   
     556  // order of magnitude of the centroid speed   
    480557
    481558  // FIXME(Ole): Work in progress
    … …  
    486563 
    487564  printf("alpha = %f, uh_i=%f, uh=%f, h_i=%f, h=%f\n",
    488         alpha, uh_i, uh, h_i, h);
     565    alpha, uh_i, uh, h_i, h);
    489566     
    490567 
    … …  
    500577  if (h < epsilon) {
    501578    b = 1.0e10; // Limit
    502   } else {     
     579  } else { 
    503580    b = m*fabs(h_i - h)/h;
    504581  }
    … …  
    507584
    508585  if (a > b) {
    509     estimate = (m-1)/(a-b);             
     586    estimate = (m-1)/(a-b);     
    510587   
    511588    printf("Alpha %f, estimate %f\n",
    512            alpha, estimate);   
    513            
     589       alpha, estimate);   
     590       
    514591    if (alpha < estimate) {
    515592      printf("Adjusting alpha from %f to %f\n",
    516              alpha, estimate);
     593         alpha, estimate);
    517594      alpha = estimate;
    518595    }
    … …  
    520597 
    521598    if (h < h_i) {
    522       estimate = (m-1)/(a-b);                
     599      estimate = (m-1)/(a-b);            
    523600   
    524601      printf("Alpha %f, estimate %f\n",
    525              alpha, estimate);   
    526            
     602         alpha, estimate);   
     603       
    527604      if (alpha < estimate) {
    528         printf("Adjusting alpha from %f to %f\n",
    529                alpha, estimate);
    530         alpha = estimate;
     605    printf("Adjusting alpha from %f to %f\n",
     606           alpha, estimate);
     607    alpha = estimate;
    531608      }   
    532609    }
    … …  
    540617
    541618int _balance_deep_and_shallow(int N,
    542                               double* wc,
    543                               double* zc,
    544                               double* wv,
    545                               double* zv,
    546                               double* hvbar, // Retire this
    547                               double* xmomc,
    548                               double* ymomc,
    549                               double* xmomv,
    550                               double* ymomv,
    551                               double H0,
    552                               int tight_slope_limiters,
    553                               int use_centroid_velocities,
    554                               double alpha_balance) {
     619                  double* wc,
     620                  double* zc,
     621                  double* wv,
     622                  double* zv,
     623                  double* hvbar, // Retire this
     624                  double* xmomc,
     625                  double* ymomc,
     626                  double* xmomv,
     627                  double* ymomv,
     628                  double H0,
     629                  int tight_slope_limiters,
     630                  int use_centroid_velocities,
     631                  double alpha_balance) {
    555632
    556633  int k, k3, i;
    … …  
    579656      // FIXME: Try with this one precomputed
    580657      for (i=0; i<3; i++) {
    581         dz = max(dz, fabs(zv[k3+i]-zc[k]));
     658    dz = max(dz, fabs(zv[k3+i]-zc[k]));
    582659      }
    583660    }
    … …  
    592669
    593670    //if (hmin < 0.0 ) {
    594     //  printf("hmin = %f\n",hmin);
     671    //  printf("hmin = %f\n",hmin);
    595672    //}
    596673   
    … …  
    607684     
    608685      if (dz > 0.0) {
    609         alpha = max( min( alpha_balance*hmin/dz, 1.0), 0.0 );     
     686    alpha = max( min( alpha_balance*hmin/dz, 1.0), 0.0 );     
    610687      } else {
    611         alpha = 1.0;  // Flat bed
     688    alpha = 1.0;  // Flat bed
    612689      }
    613690      //printf("Using old style limiter\n");
    … …  
    622699   
    623700      if (hmin < H0) {
    624         alpha = 1.0;
    625         for (i=0; i<3; i++) {
    626          
    627           h_diff = hc_k - hv[i];         
    628           if (h_diff <= 0) {
    629             // Deep water triangle is further away from bed than
    630             // shallow water (hbar < h). Any alpha will do
    631          
    632           } else { 
    633             // Denominator is positive which means that we need some of the
    634             // h-limited stage.
    635            
    636             alpha = min(alpha, (hc_k - H0)/h_diff);        
    637           }
    638         }
    639 
    640         // Ensure alpha in [0,1]
    641         if (alpha>1.0) alpha=1.0;
    642         if (alpha<0.0) alpha=0.0;
    643        
     701    alpha = 1.0;
     702    for (i=0; i<3; i++) {
     703     
     704      h_diff = hc_k - hv[i];     
     705      if (h_diff <= 0) {
     706        // Deep water triangle is further away from bed than
     707        // shallow water (hbar < h). Any alpha will do
     708     
     709      } else { 
     710        // Denominator is positive which means that we need some of the
     711        // h-limited stage.
     712       
     713        alpha = min(alpha, (hc_k - H0)/h_diff);    
     714      }
     715    }
     716
     717    // Ensure alpha in [0,1]
     718    if (alpha>1.0) alpha=1.0;
     719    if (alpha<0.0) alpha=0.0;
     720   
    644721      } else {
    645         // Use w-limited stage exclusively in deeper water.
    646         alpha = 1.0;       
     722    // Use w-limited stage exclusively in deeper water.
     723    alpha = 1.0;       
    647724      }
    648725    }
    649726
    650727
    651     //  Let
     728    //  Let
    652729    //
    653     //    wvi be the w-limited stage (wvi = zvi + hvi)
    654     //    wvi- be the h-limited state (wvi- = zvi + hvi-)
     730    //    wvi be the w-limited stage (wvi = zvi + hvi)
     731    //    wvi- be the h-limited state (wvi- = zvi + hvi-)
    655732    //
    656733    //
    657     //  where i=0,1,2 denotes the vertex ids
     734    //  where i=0,1,2 denotes the vertex ids
    658735    //
    659736    //  Weighted balance between w-limited and h-limited stage is
    660737    //
    661     //    wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi)
     738    //    wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi)
    662739    //
    663740    //  It follows that the updated wvi is
    664741    //    wvi := zvi + (1-alpha)*hvi- + alpha*hvi
    665742    //
    666     //  Momentum is balanced between constant and limited
     743    //  Momentum is balanced between constant and limited
    667744
    668745
    … …  
    670747      for (i=0; i<3; i++) {
    671748     
    672         wv[k3+i] = zv[k3+i] + (1-alpha)*hc_k + alpha*hv[i];     
    673 
    674         // Update momentum at vertices
    675         if (use_centroid_velocities == 1) {
    676           // This is a simple, efficient and robust option
    677           // It uses first order approximation of velocities, but retains
    678           // the order used by stage.
    679        
    680           // Speeds at centroids
    681           if (hc_k > epsilon) {
    682             uc = xmomc[k]/hc_k;
    683             vc = ymomc[k]/hc_k;
    684           } else {
    685             uc = 0.0;
    686             vc = 0.0;
    687           }
    688          
    689           // Vertex momenta guaranteed to be consistent with depth guaranteeing
    690           // controlled speed
    691           hv[i] = wv[k3+i] - zv[k3+i]; // Recompute (balanced) vertex depth
    692           xmomv[k3+i] = uc*hv[i];
    693           ymomv[k3+i] = vc*hv[i];
    694          
    695         } else {
    696           // Update momentum as a linear combination of
    697           // xmomc and ymomc (shallow) and momentum
    698           // from extrapolator xmomv and ymomv (deep).
    699           // This assumes that values from xmomv and ymomv have
    700           // been established e.g. by the gradient limiter.
    701 
    702           // FIXME (Ole): I think this should be used with vertex momenta
    703           // computed above using centroid_velocities instead of xmomc
    704           // and ymomc as they'll be more representative first order
    705           // values.
    706          
    707           xmomv[k3+i] = (1-alpha)*xmomc[k] + alpha*xmomv[k3+i];
    708           ymomv[k3+i] = (1-alpha)*ymomc[k] + alpha*ymomv[k3+i];
    709        
    710         }
     749    wv[k3+i] = zv[k3+i] + (1-alpha)*hc_k + alpha*hv[i];
     750
     751    // Update momentum at vertices
     752    if (use_centroid_velocities == 1) {
     753      // This is a simple, efficient and robust option
     754      // It uses first order approximation of velocities, but retains
     755      // the order used by stage.
     756   
     757      // Speeds at centroids
     758      if (hc_k > epsilon) {
     759        uc = xmomc[k]/hc_k;
     760        vc = ymomc[k]/hc_k;
     761      } else {
     762        uc = 0.0;
     763        vc = 0.0;
     764      }
     765     
     766      // Vertex momenta guaranteed to be consistent with depth guaranteeing
     767      // controlled speed
     768      hv[i] = wv[k3+i] - zv[k3+i]; // Recompute (balanced) vertex depth
     769      xmomv[k3+i] = uc*hv[i];
     770      ymomv[k3+i] = vc*hv[i];
     771     
     772    } else {
     773      // Update momentum as a linear combination of
     774      // xmomc and ymomc (shallow) and momentum
     775      // from extrapolator xmomv and ymomv (deep).
     776      // This assumes that values from xmomv and ymomv have
     777      // been established e.g. by the gradient limiter.
     778
     779      // FIXME (Ole): I think this should be used with vertex momenta
     780      // computed above using centroid_velocities instead of xmomc
     781      // and ymomc as they'll be more representative first order
     782      // values.
     783     
     784      xmomv[k3+i] = (1-alpha)*xmomc[k] + alpha*xmomv[k3+i];
     785      ymomv[k3+i] = (1-alpha)*ymomc[k] + alpha*ymomv[k3+i];
     786   
     787    }
    711788      }
    712789    }
    … …  
    719796
    720797int _protect(int N,
    721              double minimum_allowed_height,
    722              double maximum_allowed_speed,
    723              double epsilon,
    724              double* wc,
    725              double* zc,
    726              double* xmomc,
    727              double* ymomc) {
     798         double minimum_allowed_height,
     799         double maximum_allowed_speed,
     800         double epsilon,
     801         double* wc,
     802         double* zc,
     803         double* xmomc,
     804         double* ymomc) {
    728805
    729806  int k;
    … …  
    738815
    739816      if (hc < minimum_allowed_height) {
    740        
    741         // Set momentum to zero and ensure h is non negative
    742         xmomc[k] = 0.0;
    743         ymomc[k] = 0.0;
    744         if (hc <= 0.0) wc[k] = zc[k];
     817       
     818    // Set momentum to zero and ensure h is non negative
     819    xmomc[k] = 0.0;
     820    ymomc[k] = 0.0;
     821    if (hc <= 0.0) wc[k] = zc[k];
    745822      }
    746823    }
    … …  
    757834             
    758835        if (hc <= 0.0) {
    759                 wc[k] = zc[k];
    760         xmomc[k] = 0.0;
    761         ymomc[k] = 0.0;
     836            wc[k] = zc[k];
     837        xmomc[k] = 0.0;
     838        ymomc[k] = 0.0;
    762839        } else {
    763840          //Reduce excessive speeds derived from division by small hc
    764         //FIXME (Ole): This may be unnecessary with new slope limiters
    765         //in effect.
     841        //FIXME (Ole): This may be unnecessary with new slope limiters
     842        //in effect.
    766843         
    767844          u = xmomc[k]/hc;
    768           if (fabs(u) > maximum_allowed_speed) {
    769             reduced_speed = maximum_allowed_speed * u/fabs(u);
    770             //printf("Speed (u) has been reduced from %.3f to %.3f\n",
    771             //  u, reduced_speed);
    772             xmomc[k] = reduced_speed * hc;
    773           }
     845      if (fabs(u) > maximum_allowed_speed) {
     846        reduced_speed = maximum_allowed_speed * u/fabs(u);
     847        //printf("Speed (u) has been reduced from %.3f to %.3f\n",
     848        //  u, reduced_speed);
     849        xmomc[k] = reduced_speed * hc;
     850      }
    774851
    775852          v = ymomc[k]/hc;
    776           if (fabs(v) > maximum_allowed_speed) {
    777             reduced_speed = maximum_allowed_speed * v/fabs(v);
    778             //printf("Speed (v) has been reduced from %.3f to %.3f\n",
    779             //  v, reduced_speed);
    780             ymomc[k] = reduced_speed * hc;
    781           }
     853      if (fabs(v) > maximum_allowed_speed) {
     854        reduced_speed = maximum_allowed_speed * v/fabs(v);
     855        //printf("Speed (v) has been reduced from %.3f to %.3f\n",
     856        //  v, reduced_speed);
     857        ymomc[k] = reduced_speed * hc;
     858      }
    782859        }
    783860      }
    … …  
    790867
    791868int _assign_wind_field_values(int N,
    792                               double* xmom_update,
    793                               double* ymom_update,
    794                               double* s_vec,
    795                               double* phi_vec,
    796                               double cw) {
     869                  double* xmom_update,
     870                  double* ymom_update,
     871                  double* s_vec,
     872                  double* phi_vec,
     873                  double cw) {
    797874
    798875  // Assign windfield values to momentum updates
    … …  
    839916
    840917  if (!PyArg_ParseTuple(args, "OOOddOddd",
    841                         &normal, &ql, &qr, &zl, &zr, &edgeflux,
    842                         &epsilon, &g, &H0)) {
     918            &normal, &ql, &qr, &zl, &zr, &edgeflux,
     919            &epsilon, &g, &H0)) {
    843920    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: flux_function_central could not parse input arguments");
    844921    return NULL;
    … …  
    847924 
    848925  _flux_function_central((double*) ql -> data,
    849                         (double*) qr -> data,
    850                         zl,
    851                          zr,                                           
    852                         ((double*) normal -> data)[0],
    853                          ((double*) normal -> data)[1],                 
    854                         epsilon, H0, g,
    855                         (double*) edgeflux -> data,
    856                         &max_speed);
     926            (double*) qr -> data,
     927            zl,
     928             zr,                       
     929            ((double*) normal -> data)[0],
     930             ((double*) normal -> data)[1],         
     931            epsilon, H0, g,
     932            (double*) edgeflux -> data,
     933            &max_speed);
    857934 
    858935  return Py_BuildValue("d", max_speed); 
    … …  
    874951
    875952  if (!PyArg_ParseTuple(args, "dOOOOO",
    876                         &g, &h, &v, &x,
    877                         &xmom, &ymom)) {
     953            &g, &h, &v, &x,
     954            &xmom, &ymom)) {
    878955    //&epsilon)) {
    879956    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: gravity could not parse input arguments");
    … …  
    9401017
    9411018  if (!PyArg_ParseTuple(args, "ddOOOOOOO",
    942                         &g, &eps, &w, &z, &uh, &vh, &eta,
    943                         &xmom, &ymom)) {
     1019            &g, &eps, &w, &z, &uh, &vh, &eta,
     1020            &xmom, &ymom)) {
    9441021    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: manning_friction could not parse input arguments");
    9451022    return NULL;
    … …  
    9571034  N = w -> dimensions[0];
    9581035  _manning_friction(g, eps, N,
    959                     (double*) w -> data,
    960                     (double*) z -> data,
    961                     (double*) uh -> data,
    962                     (double*) vh -> data,
    963                     (double*) eta -> data,
    964                     (double*) xmom -> data,
    965                     (double*) ymom -> data);
     1036            (double*) w -> data,
     1037            (double*) z -> data,
     1038            (double*) uh -> data,
     1039            (double*) vh -> data,
     1040            (double*) eta -> data,
     1041            (double*) xmom -> data,
     1042            (double*) ymom -> data);
    9661043
    9671044  return Py_BuildValue("");
    … …  
    9811058
    9821059  if (!PyArg_ParseTuple(args, "ddOOOOOOO",
    983                         &g, &eps, &w, &z, &uh, &vh, &eta,
    984                         &xmom, &ymom))
     1060            &g, &eps, &w, &z, &uh, &vh, &eta,
     1061            &xmom, &ymom))
    9851062    return NULL;
    9861063
    … …  
    9961073  N = w -> dimensions[0];
    9971074  _manning_friction_explicit(g, eps, N,
    998                     (double*) w -> data,
    999                     (double*) z -> data,
    1000                     (double*) uh -> data,
    1001                     (double*) vh -> data,
    1002                     (double*) eta -> data,
    1003                     (double*) xmom -> data,
    1004                     (double*) ymom -> data);
     1075            (double*) w -> data,
     1076            (double*) z -> data,
     1077            (double*) uh -> data,
     1078            (double*) vh -> data,
     1079            (double*) eta -> data,
     1080            (double*) xmom -> data,
     1081            (double*) ymom -> data);
    10051082
    10061083  return Py_BuildValue("");
    … …  
    10121089// Computational routine
    10131090int _extrapolate_second_order_sw(int number_of_elements,
    1014                                   double epsilon,
    1015                                   double minimum_allowed_height,
    1016                                   double beta_w,
    1017                                   double beta_w_dry,
    1018                                   double beta_uh,
    1019                                   double beta_uh_dry,
    1020                                   double beta_vh,
    1021                                   double beta_vh_dry,
    1022                                   long* surrogate_neighbours,
    1023                                   long* number_of_boundaries,
    1024                                   double* centroid_coordinates,
    1025                                   double* stage_centroid_values,
    1026                                   double* xmom_centroid_values,
    1027                                   double* ymom_centroid_values,
    1028                                   double* elevation_centroid_values,
    1029                                   double* vertex_coordinates,
    1030                                   double* stage_vertex_values,
    1031                                   double* xmom_vertex_values,
    1032                                   double* ymom_vertex_values,
    1033                                   double* elevation_vertex_values,
    1034                                   int optimise_dry_cells) {
    1035                                  
    1036                                  
     1091                                 double epsilon,
     1092                                 double minimum_allowed_height,
     1093                                 double beta_w,
     1094                                 double beta_w_dry,
     1095                                 double beta_uh,
     1096                                 double beta_uh_dry,
     1097                                 double beta_vh,
     1098                                 double beta_vh_dry,
     1099                                 long* surrogate_neighbours,
     1100                                 long* number_of_boundaries,
     1101                                 double* centroid_coordinates,
     1102                                 double* stage_centroid_values,
     1103                                 double* xmom_centroid_values,
     1104                                 double* ymom_centroid_values,
     1105                                 double* elevation_centroid_values,
     1106                                 double* vertex_coordinates,
     1107                                 double* stage_vertex_values,
     1108                                 double* xmom_vertex_values,
     1109                                 double* ymom_vertex_values,
     1110                                 double* elevation_vertex_values,
     1111                                 int optimise_dry_cells) {
     1112                 
     1113                 
    10371114
    10381115  // Local variables
    10391116  double a, b; // Gradient vector used to calculate vertex values from centroids
    1040   int k,k0,k1,k2,k3,k6,coord_index,i;
    1041   double x,y,x0,y0,x1,y1,x2,y2,xv0,yv0,xv1,yv1,xv2,yv2; // Vertices of the auxiliary triangle
    1042   double dx1,dx2,dy1,dy2,dxv0,dxv1,dxv2,dyv0,dyv1,dyv2,dq0,dq1,dq2,area2;
     1117  int k, k0, k1, k2, k3, k6, coord_index, i;
     1118  double x, y, x0, y0, x1, y1, x2, y2, xv0, yv0, xv1, yv1, xv2, yv2; // Vertices of the auxiliary triangle
     1119  double dx1, dx2, dy1, dy2, dxv0, dxv1, dxv2, dyv0, dyv1, dyv2, dq0, dq1, dq2, area2, inv_area2;
    10431120  double dqv[3], qmin, qmax, hmin, hmax;
    10441121  double hc, h0, h1, h2, beta_tmp, hfactor;
    10451122
    1046        
    1047   for (k=0; k<number_of_elements; k++) {
     1123   
     1124  for (k = 0; k < number_of_elements; k++)
     1125  {
    10481126    k3=k*3;
    10491127    k6=k*6;
    10501128
    1051    
    1052     if (number_of_boundaries[k]==3){
     1129    if (number_of_boundaries[k]==3)
     1130    {
    10531131      // No neighbours, set gradient on the triangle to zero
    10541132     
    … …  
    10651143      continue;
    10661144    }
    1067     else {
     1145    else
     1146    {
    10681147      // Triangle k has one or more neighbours.
    10691148      // Get centroid and vertex coordinates of the triangle
    10701149     
    10711150      // Get the vertex coordinates
    1072       xv0 = vertex_coordinates[k6];   yv0=vertex_coordinates[k6+1];
    1073       xv1 = vertex_coordinates[k6+2]; yv1=vertex_coordinates[k6+3];
    1074       xv2 = vertex_coordinates[k6+4]; yv2=vertex_coordinates[k6+5];
     1151      xv0 = vertex_coordinates[k6];   
     1152      yv0 = vertex_coordinates[k6+1];
     1153      xv1 = vertex_coordinates[k6+2];
     1154      yv1 = vertex_coordinates[k6+3];
     1155      xv2 = vertex_coordinates[k6+4];
     1156      yv2 = vertex_coordinates[k6+5];
    10751157     
    10761158      // Get the centroid coordinates
    1077       coord_index=2*k;
    1078       x=centroid_coordinates[coord_index];
    1079       y=centroid_coordinates[coord_index+1];
     1159      coord_index = 2*k;
     1160      x = centroid_coordinates[coord_index];
     1161      y = centroid_coordinates[coord_index+1];
    10801162     
    10811163      // Store x- and y- differentials for the vertices of
    10821164      // triangle k relative to the centroid
    1083       dxv0=xv0-x; dxv1=xv1-x; dxv2=xv2-x;
    1084       dyv0=yv0-y; dyv1=yv1-y; dyv2=yv2-y;
     1165      dxv0 = xv0 - x;
     1166      dxv1 = xv1 - x;
     1167      dxv2 = xv2 - x;
     1168      dyv0 = yv0 - y;
     1169      dyv1 = yv1 - y;
     1170      dyv2 = yv2 - y;
    10851171    }
    10861172
    10871173
    10881174           
    1089     if (number_of_boundaries[k]<=1){
    1090    
     1175    if (number_of_boundaries[k]<=1)
     1176    {
    10911177      //==============================================
    10921178      // Number of boundaries <= 1
    … …  
    10991185      // from this centroid and its two neighbours
    11001186     
    1101       k0=surrogate_neighbours[k3];
    1102       k1=surrogate_neighbours[k3+1];
    1103       k2=surrogate_neighbours[k3+2];
     1187      k0 = surrogate_neighbours[k3];
     1188      k1 = surrogate_neighbours[k3 + 1];
     1189      k2 = surrogate_neighbours[k3 + 2];
    11041190     
    11051191      // Get the auxiliary triangle's vertex coordinates
    11061192      // (really the centroids of neighbouring triangles)
    1107       coord_index=2*k0;
    1108       x0=centroid_coordinates[coord_index];
    1109       y0=centroid_coordinates[coord_index+1];
    1110      
    1111       coord_index=2*k1;
    1112       x1=centroid_coordinates[coord_index];
    1113       y1=centroid_coordinates[coord_index+1];
    1114      
    1115       coord_index=2*k2;
    1116       x2=centroid_coordinates[coord_index];
    1117       y2=centroid_coordinates[coord_index+1];
     1193      coord_index = 2*k0;
     1194      x0 = centroid_coordinates[coord_index];
     1195      y0 = centroid_coordinates[coord_index+1];
     1196     
     1197      coord_index = 2*k1;
     1198      x1 = centroid_coordinates[coord_index];
     1199      y1 = centroid_coordinates[coord_index+1];
     1200     
     1201      coord_index = 2*k2;
     1202      x2 = centroid_coordinates[coord_index];
     1203      y2 = centroid_coordinates[coord_index+1];
    11181204     
    11191205      // Store x- and y- differentials for the vertices
    11201206      // of the auxiliary triangle
    1121       dx1=x1-x0; dx2=x2-x0;
    1122       dy1=y1-y0; dy2=y2-y0;
     1207      dx1 = x1 - x0;
     1208      dx2 = x2 - x0;
     1209      dy1 = y1 - y0;
     1210      dy2 = y2 - y0;
    11231211     
    11241212      // Calculate 2*area of the auxiliary triangle
    … …  
    11281216      // If the mesh is 'weird' near the boundary,
    11291217      // the triangle might be flat or clockwise:
    1130       if (area2<=0) {
    1131         PyErr_SetString(PyExc_RuntimeError,
    1132                         "shallow_water_ext.c: negative triangle area encountered");
    1133         return -1;
     1218      if (area2 <= 0)
     1219      {
     1220        PyErr_SetString(PyExc_RuntimeError,
     1221                        "shallow_water_ext.c: negative triangle area encountered");
     1222   
     1223        return -1;
    11341224      } 
    11351225     
    … …  
    11391229      h1 = stage_centroid_values[k1] - elevation_centroid_values[k1];
    11401230      h2 = stage_centroid_values[k2] - elevation_centroid_values[k2];
    1141       hmin = min(min(h0,min(h1,h2)),hc);
     1231      hmin = min(min(h0, min(h1, h2)), hc);
    11421232      //hfactor = hc/(hc + 1.0);
    11431233
    11441234      hfactor = 0.0;
    1145       if (hmin > 0.001 ) {
    1146           hfactor = (hmin-0.001)/(hmin+0.004);
    1147       }
    1148      
    1149       if (optimise_dry_cells) {     
    1150         // Check if linear reconstruction is necessary for triangle k
    1151         // This check will exclude dry cells.
    1152 
    1153         hmax = max(h0,max(h1,h2));     
    1154         if (hmax < epsilon) {
    1155           continue;
    1156         }
    1157       }
    1158 
    1159            
     1235      if (hmin > 0.001)
     1236      {
     1237        hfactor = (hmin - 0.001)/(hmin + 0.004);
     1238      }
     1239     
     1240      if (optimise_dry_cells)
     1241      {     
     1242        // Check if linear reconstruction is necessary for triangle k
     1243        // This check will exclude dry cells.
     1244
     1245        hmax = max(h0, max(h1, h2));     
     1246        if (hmax < epsilon)
     1247        {
     1248            continue;
     1249        }
     1250      }
     1251   
    11601252      //-----------------------------------
    11611253      // stage
    … …  
    11641256      // Calculate the difference between vertex 0 of the auxiliary
    11651257      // triangle and the centroid of triangle k
    1166       dq0=stage_centroid_values[k0]-stage_centroid_values[k];
     1258      dq0 = stage_centroid_values[k0] - stage_centroid_values[k];
    11671259     
    11681260      // Calculate differentials between the vertices
    11691261      // of the auxiliary triangle (centroids of neighbouring triangles)
    1170       dq1=stage_centroid_values[k1]-stage_centroid_values[k0];
    1171       dq2=stage_centroid_values[k2]-stage_centroid_values[k0];
    1172      
     1262      dq1 = stage_centroid_values[k1] - stage_centroid_values[k0];
     1263      dq2 = stage_centroid_values[k2] - stage_centroid_values[k0];
     1264     
     1265          inv_area2 = 1.0/area2;
    11731266      // Calculate the gradient of stage on the auxiliary triangle
    11741267      a = dy2*dq1 - dy1*dq2;
    1175       a /= area2;
     1268      a *= inv_area2;
    11761269      b = dx1*dq2 - dx2*dq1;
    1177       b /= area2;
     1270      b *= inv_area2;
    11781271     
    11791272      // Calculate provisional jumps in stage from the centroid
    11801273      // of triangle k to its vertices, to be limited
    1181       dqv[0]=a*dxv0+b*dyv0;
    1182       dqv[1]=a*dxv1+b*dyv1;
    1183       dqv[2]=a*dxv2+b*dyv2;
     1274      dqv[0] = a*dxv0 + b*dyv0;
     1275      dqv[1] = a*dxv1 + b*dyv1;
     1276      dqv[2] = a*dxv2 + b*dyv2;
    11841277     
    11851278      // Now we want to find min and max of the centroid and the
    11861279      // vertices of the auxiliary triangle and compute jumps
    11871280      // from the centroid to the min and max
    1188       find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
     1281      find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
    11891282     
    11901283      // Playing with dry wet interface
    … …  
    11931286      //if (hmin>minimum_allowed_height)
    11941287      beta_tmp = beta_w_dry + (beta_w - beta_w_dry) * hfactor;
    1195        
     1288   
    11961289      //printf("min_alled_height = %f\n",minimum_allowed_height);
    11971290      //printf("hmin = %f\n",hmin);
    … …  
    11991292      //printf("beta_tmp = %f\n",beta_tmp);
    12001293      // Limit the gradient
    1201       limit_gradient(dqv,qmin,qmax,beta_tmp);
    1202      
    1203       for (i=0;i<3;i++)
    1204         stage_vertex_values[k3+i]=stage_centroid_values[k]+dqv[i];
     1294      limit_gradient(dqv, qmin, qmax, beta_tmp);
     1295     
     1296      //for (i=0;i<3;i++)
     1297      stage_vertex_values[k3+0] = stage_centroid_values[k] + dqv[0];
     1298          stage_vertex_values[k3+1] = stage_centroid_values[k] + dqv[1];
     1299          stage_vertex_values[k3+2] = stage_centroid_values[k] + dqv[2];
    12051300     
    12061301     
    … …  
    12111306      // Calculate the difference between vertex 0 of the auxiliary
    12121307      // triangle and the centroid of triangle k     
    1213       dq0=xmom_centroid_values[k0]-xmom_centroid_values[k];
     1308      dq0 = xmom_centroid_values[k0] - xmom_centroid_values[k];
    12141309     
    12151310      // Calculate differentials between the vertices
    12161311      // of the auxiliary triangle
    1217       dq1=xmom_centroid_values[k1]-xmom_centroid_values[k0];
    1218       dq2=xmom_centroid_values[k2]-xmom_centroid_values[k0];
     1312      dq1 = xmom_centroid_values[k1] - xmom_centroid_values[k0];
     1313      dq2 = xmom_centroid_values[k2] - xmom_centroid_values[k0];
    12191314     
    12201315      // Calculate the gradient of xmom on the auxiliary triangle
    12211316      a = dy2*dq1 - dy1*dq2;
    1222       a /= area2;
     1317      a *= inv_area2;
    12231318      b = dx1*dq2 - dx2*dq1;
    1224       b /= area2;
     1319      b *= inv_area2;
    12251320     
    12261321      // Calculate provisional jumps in stage from the centroid
    12271322      // of triangle k to its vertices, to be limited     
    1228       dqv[0]=a*dxv0+b*dyv0;
    1229       dqv[1]=a*dxv1+b*dyv1;
    1230       dqv[2]=a*dxv2+b*dyv2;
     1323      dqv[0] = a*dxv0+b*dyv0;
     1324      dqv[1] = a*dxv1+b*dyv1;
     1325      dqv[2] = a*dxv2+b*dyv2;
    12311326     
    12321327      // Now we want to find min and max of the centroid and the
    12331328      // vertices of the auxiliary triangle and compute jumps
    12341329      // from the centroid to the min and max
    1235       find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
     1330      find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
    12361331      //beta_tmp = beta_uh;
    12371332      //if (hmin<minimum_allowed_height)
    … …  
    12401335
    12411336      // Limit the gradient
    1242       limit_gradient(dqv,qmin,qmax,beta_tmp);
    1243 
    1244       for (i=0;i<3;i++)
    1245         xmom_vertex_values[k3+i]=xmom_centroid_values[k]+dqv[i];
    1246      
     1337      limit_gradient(dqv, qmin, qmax, beta_tmp);
     1338
     1339      for (i=0; i < 3; i++)
     1340      {
     1341        xmom_vertex_values[k3+i] = xmom_centroid_values[k] + dqv[i];
     1342      }
    12471343     
    12481344      //-----------------------------------
    … …  
    12521348      // Calculate the difference between vertex 0 of the auxiliary
    12531349      // triangle and the centroid of triangle k
    1254       dq0=ymom_centroid_values[k0]-ymom_centroid_values[k];
     1350      dq0 = ymom_centroid_values[k0] - ymom_centroid_values[k];
    12551351     
    12561352      // Calculate differentials between the vertices
    12571353      // of the auxiliary triangle
    1258       dq1=ymom_centroid_values[k1]-ymom_centroid_values[k0];
    1259       dq2=ymom_centroid_values[k2]-ymom_centroid_values[k0];
     1354      dq1 = ymom_centroid_values[k1] - ymom_centroid_values[k0];
     1355      dq2 = ymom_centroid_values[k2] - ymom_centroid_values[k0];
    12601356     
    12611357      // Calculate the gradient of xmom on the auxiliary triangle
    12621358      a = dy2*dq1 - dy1*dq2;
    1263       a /= area2;
     1359      a *= inv_area2;
    12641360      b = dx1*dq2 - dx2*dq1;
    1265       b /= area2;
     1361      b *= inv_area2;
    12661362     
    12671363      // Calculate provisional jumps in stage from the centroid
    12681364      // of triangle k to its vertices, to be limited
    1269       dqv[0]=a*dxv0+b*dyv0;
    1270       dqv[1]=a*dxv1+b*dyv1;
    1271       dqv[2]=a*dxv2+b*dyv2;
     1365      dqv[0] = a*dxv0 + b*dyv0;
     1366      dqv[1] = a*dxv1 + b*dyv1;
     1367      dqv[2] = a*dxv2 + b*dyv2;
    12721368     
    12731369      // Now we want to find min and max of the centroid and the
    12741370      // vertices of the auxiliary triangle and compute jumps
    12751371      // from the centroid to the min and max
    1276       find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
     1372      find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
    12771373     
    12781374      //beta_tmp = beta_vh;
    … …  
    12801376      //if (hmin<minimum_allowed_height)
    12811377      //beta_tmp = beta_vh_dry;
    1282       beta_tmp = beta_vh_dry + (beta_vh - beta_vh_dry) * hfactor;       
     1378      beta_tmp = beta_vh_dry + (beta_vh - beta_vh_dry) * hfactor;   
    12831379
    12841380      // Limit the gradient
    1285       limit_gradient(dqv,qmin,qmax,beta_tmp);
     1381      limit_gradient(dqv, qmin, qmax, beta_tmp);
    12861382     
    12871383      for (i=0;i<3;i++)
    1288         ymom_vertex_values[k3+i]=ymom_centroid_values[k]+dqv[i];
    1289        
     1384      {
     1385        ymom_vertex_values[k3 + i] = ymom_centroid_values[k] + dqv[i];
     1386      }
    12901387    } // End number_of_boundaries <=1
    1291     else{
     1388    else
     1389    {
    12921390
    12931391      //==============================================
    … …  
    12981396     
    12991397      // Find the only internal neighbour (k1?)
    1300       for (k2=k3;k2<k3+3;k2++){
    1301         // Find internal neighbour of triangle k     
    1302         // k2 indexes the edges of triangle k   
    1303      
    1304         if (surrogate_neighbours[k2]!=k)
    1305           break;
    1306       }
    1307      
    1308       if ((k2==k3+3)) {
    1309         // If we didn't find an internal neighbour
    1310         PyErr_SetString(PyExc_RuntimeError,
    1311                         "shallow_water_ext.c: Internal neighbour not found");     
    1312         return -1;
    1313       }
    1314      
    1315       k1=surrogate_neighbours[k2];
     1398      for (k2 = k3; k2 < k3 + 3; k2++)
     1399      {
     1400      // Find internal neighbour of triangle k     
     1401      // k2 indexes the edges of triangle k   
     1402     
     1403          if (surrogate_neighbours[k2] != k)
     1404          {
     1405             break;
     1406          }
     1407      }
     1408     
     1409      if ((k2 == k3 + 3))
     1410      {
     1411        // If we didn't find an internal neighbour
     1412        PyErr_SetString(PyExc_RuntimeError,
     1413                        "shallow_water_ext.c: Internal neighbour not found");     
     1414        return -1;
     1415      }
     1416     
     1417      k1 = surrogate_neighbours[k2];
    13161418     
    13171419      // The coordinates of the triangle are already (x,y).
    13181420      // Get centroid of the neighbour (x1,y1)
    1319       coord_index=2*k1;
    1320       x1=centroid_coordinates[coord_index];
    1321       y1=centroid_coordinates[coord_index+1];
     1421      coord_index = 2*k1;
     1422      x1 = centroid_coordinates[coord_index];
     1423      y1 = centroid_coordinates[coord_index + 1];
    13221424     
    13231425      // Compute x- and y- distances between the centroid of
    13241426      // triangle k and that of its neighbour
    1325       dx1=x1-x; dy1=y1-y;
     1427      dx1 = x1 - x;
     1428      dy1 = y1 - y;
    13261429     
    13271430      // Set area2 as the square of the distance
    1328       area2=dx1*dx1+dy1*dy1;
     1431      area2 = dx1*dx1 + dy1*dy1;
    13291432     
    13301433      // Set dx2=(x1-x0)/((x1-x0)^2+(y1-y0)^2)
    … …  
    13321435      // respectively correspond to the x- and y- gradients
    13331436      // of the conserved quantities
    1334       dx2=1.0/area2;
    1335       dy2=dx2*dy1;
    1336       dx2*=dx1;
     1437      dx2 = 1.0/area2;
     1438      dy2 = dx2*dy1;
     1439      dx2 *= dx1;
    13371440     
    13381441     
    … …  
    13421445
    13431446      // Compute differentials
    1344       dq1=stage_centroid_values[k1]-stage_centroid_values[k];
     1447      dq1 = stage_centroid_values[k1] - stage_centroid_values[k];
    13451448     
    13461449      // Calculate the gradient between the centroid of triangle k
    13471450      // and that of its neighbour
    1348       a=dq1*dx2;
    1349       b=dq1*dy2;
     1451      a = dq1*dx2;
     1452      b = dq1*dy2;
    13501453     
    13511454      // Calculate provisional vertex jumps, to be limited
    1352       dqv[0]=a*dxv0+b*dyv0;
    1353       dqv[1]=a*dxv1+b*dyv1;
    1354       dqv[2]=a*dxv2+b*dyv2;
     1455      dqv[0] = a*dxv0 + b*dyv0;
     1456      dqv[1] = a*dxv1 + b*dyv1;
     1457      dqv[2] = a*dxv2 + b*dyv2;
    13551458     
    13561459      // Now limit the jumps
    1357       if (dq1>=0.0){
    1358         qmin=0.0;
    1359         qmax=dq1;
    1360       }
    1361       else{
    1362         qmin=dq1;
    1363         qmax=0.0;
     1460      if (dq1>=0.0)
     1461      {
     1462        qmin=0.0;
     1463        qmax=dq1;
     1464      }
     1465      else
     1466      {
     1467        qmin = dq1;
     1468        qmax = 0.0;
    13641469      }
    13651470     
    13661471      // Limit the gradient
    1367       limit_gradient(dqv,qmin,qmax,beta_w);
    1368      
    1369       for (i=0;i<3;i++)
    1370         stage_vertex_values[k3+i]=stage_centroid_values[k]+dqv[i];
    1371      
     1472      limit_gradient(dqv, qmin, qmax, beta_w);
     1473     
     1474      //for (i=0; i < 3; i++)
     1475      //{
     1476      stage_vertex_values[k3] = stage_centroid_values[k] + dqv[0];
     1477      stage_vertex_values[k3 + 1] = stage_centroid_values[k] + dqv[1];
     1478      stage_vertex_values[k3 + 2] = stage_centroid_values[k] + dqv[2];
     1479      //}
    13721480     
    13731481      //-----------------------------------
    … …  
    13761484     
    13771485      // Compute differentials
    1378       dq1=xmom_centroid_values[k1]-xmom_centroid_values[k];
     1486      dq1 = xmom_centroid_values[k1] - xmom_centroid_values[k];
    13791487     
    13801488      // Calculate the gradient between the centroid of triangle k
    13811489      // and that of its neighbour
    1382       a=dq1*dx2;
    1383       b=dq1*dy2;
     1490      a = dq1*dx2;
     1491      b = dq1*dy2;
    13841492     
    13851493      // Calculate provisional vertex jumps, to be limited
    1386       dqv[0]=a*dxv0+b*dyv0;
    1387       dqv[1]=a*dxv1+b*dyv1;
    1388       dqv[2]=a*dxv2+b*dyv2;
     1494      dqv[0] = a*dxv0+b*dyv0;
     1495      dqv[1] = a*dxv1+b*dyv1;
     1496      dqv[2] = a*dxv2+b*dyv2;
    13891497     
    13901498      // Now limit the jumps
    1391       if (dq1>=0.0){
    1392         qmin=0.0;
    1393         qmax=dq1;
    1394       }
    1395       else{
    1396         qmin=dq1;
    1397         qmax=0.0;
     1499      if (dq1 >= 0.0)
     1500      {
     1501        qmin = 0.0;
     1502        qmax = dq1;
     1503      }
     1504      else
     1505      {
     1506        qmin = dq1;
     1507        qmax = 0.0;
    13981508      }
    13991509     
    14001510      // Limit the gradient     
    1401       limit_gradient(dqv,qmin,qmax,beta_w);
    1402      
    1403       for (i=0;i<3;i++)
    1404         xmom_vertex_values[k3+i]=xmom_centroid_values[k]+dqv[i];
    1405      
     1511      limit_gradient(dqv, qmin, qmax, beta_w);
     1512     
     1513      //for (i=0;i<3;i++)
     1514      //xmom_vertex_values[k3] = xmom_centroid_values[k] + dqv[0];
     1515      //xmom_vertex_values[k3 + 1] = xmom_centroid_values[k] + dqv[1];
     1516      //xmom_vertex_values[k3 + 2] = xmom_centroid_values[k] + dqv[2];
     1517     
     1518      for (i = 0; i < 3;i++)
     1519      {
     1520          xmom_vertex_values[k3 + i] = xmom_centroid_values[k] + dqv[i];
     1521      }
    14061522     
    14071523      //-----------------------------------
    … …  
    14101526
    14111527      // Compute differentials
    1412       dq1=ymom_centroid_values[k1]-ymom_centroid_values[k];
     1528      dq1 = ymom_centroid_values[k1] - ymom_centroid_values[k];
    14131529     
    14141530      // Calculate the gradient between the centroid of triangle k
    14151531      // and that of its neighbour
    1416       a=dq1*dx2;
    1417       b=dq1*dy2;
     1532      a = dq1*dx2;
     1533      b = dq1*dy2;
    14181534     
    14191535      // Calculate provisional vertex jumps, to be limited
    1420       dqv[0]=a*dxv0+b*dyv0;
    1421       dqv[1]=a*dxv1+b*dyv1;
    1422       dqv[2]=a*dxv2+b*dyv2;
     1536      dqv[0] = a*dxv0 + b*dyv0;
     1537      dqv[1] = a*dxv1 + b*dyv1;
     1538      dqv[2] = a*dxv2 + b*dyv2;
    14231539     
    14241540      // Now limit the jumps
    1425       if (dq1>=0.0){
    1426         qmin=0.0;
    1427         qmax=dq1;
    1428       }
    1429       else{
    1430         qmin=dq1;
    1431         qmax=0.0;
     1541      if (dq1>=0.0)
     1542      {
     1543        qmin = 0.0;
     1544        qmax = dq1;
     1545      }
     1546      else
     1547      {
     1548        qmin = dq1;
     1549        qmax = 0.0;
    14321550      }
    14331551     
    14341552      // Limit the gradient
    1435       limit_gradient(dqv,qmin,qmax,beta_w);
    1436      
    1437       for (i=0;i<3;i++)
    1438         ymom_vertex_values[k3+i]=ymom_centroid_values[k]+dqv[i];
    1439        
     1553      limit_gradient(dqv, qmin, qmax, beta_w);
     1554     
     1555      //for (i=0;i<3;i++)
     1556      //ymom_vertex_values[k3] = ymom_centroid_values[k] + dqv[0];
     1557      //ymom_vertex_values[k3 + 1] = ymom_centroid_values[k] + dqv[1];
     1558      //ymom_vertex_values[k3 + 2] = ymom_centroid_values[k] + dqv[2];
     1559
     1560for (i=0;i<3;i++)
     1561        {
     1562ymom_vertex_values[k3 + i] = ymom_centroid_values[k] + dqv[i];
     1563        }
     1564//ymom_vertex_values[k3] = ymom_centroid_values[k] + dqv[0];
     1565//ymom_vertex_values[k3 + 1] = ymom_centroid_values[k] + dqv[1];
     1566//ymom_vertex_values[k3 + 2] = ymom_centroid_values[k] + dqv[2];
    14401567    } // else [number_of_boundaries==2]
    14411568  } // for k=0 to number_of_elements-1
    14421569 
    14431570  return 0;
    1444 }                       
     1571}           
    14451572
    14461573
    … …  
    14771604    Post conditions:
    14781605            The vertices of each triangle have values from a
    1479             limited linear reconstruction
    1480             based on centroid values
     1606        limited linear reconstruction
     1607        based on centroid values
    14811608
    14821609  */
    … …  
    15051632  // Convert Python arguments to C
    15061633  if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOi",
    1507                         &domain,
    1508                         &surrogate_neighbours,
    1509                         &number_of_boundaries,
    1510                         &centroid_coordinates,
    1511                         &stage_centroid_values,
    1512                         &xmom_centroid_values,
    1513                         &ymom_centroid_values,
    1514                         &elevation_centroid_values,
    1515                         &vertex_coordinates,
    1516                         &stage_vertex_values,
    1517                         &xmom_vertex_values,
    1518                         &ymom_vertex_values,
    1519                         &elevation_vertex_values,
    1520                         &optimise_dry_cells)) {                 
    1521                        
     1634            &domain,
     1635            &surrogate_neighbours,
     1636            &number_of_boundaries,
     1637            &centroid_coordinates,
     1638            &stage_centroid_values,
     1639            &xmom_centroid_values,
     1640            &ymom_centroid_values,
     1641            &elevation_centroid_values,
     1642            &vertex_coordinates,
     1643            &stage_vertex_values,
     1644            &xmom_vertex_values,
     1645            &ymom_vertex_values,
     1646            &elevation_vertex_values,
     1647            &optimise_dry_cells)) {         
     1648           
    15221649    PyErr_SetString(PyExc_RuntimeError,
    1523                     "Input arguments to extrapolate_second_order_sw failed");
     1650            "Input arguments to extrapolate_second_order_sw failed");
    15241651    return NULL;
    15251652  }
    … …  
    15571684  // Call underlying computational routine
    15581685  e = _extrapolate_second_order_sw(number_of_elements,
    1559                                    epsilon,
    1560                                    minimum_allowed_height,
    1561                                    beta_w,
    1562                                    beta_w_dry,
    1563                                    beta_uh,
    1564                                    beta_uh_dry,
    1565                                    beta_vh,
    1566                                    beta_vh_dry,
    1567                                    (long*) surrogate_neighbours -> data,
    1568                                    (long*) number_of_boundaries -> data,
    1569                                    (double*) centroid_coordinates -> data,
    1570                                    (double*) stage_centroid_values -> data,
    1571                                    (double*) xmom_centroid_values -> data,
    1572                                    (double*) ymom_centroid_values -> data,
    1573                                    (double*) elevation_centroid_values -> data,
    1574                                    (double*) vertex_coordinates -> data,
    1575                                    (double*) stage_vertex_values -> data,
    1576                                    (double*) xmom_vertex_values -> data,
    1577                                    (double*) ymom_vertex_values -> data,
    1578                                    (double*) elevation_vertex_values -> data,
    1579                                    optimise_dry_cells);
     1686                   epsilon,
     1687                   minimum_allowed_height,
     1688                   beta_w,
     1689                   beta_w_dry,
     1690                   beta_uh,
     1691                   beta_uh_dry,
     1692                   beta_vh,
     1693                   beta_vh_dry,
     1694                   (long*) surrogate_neighbours -> data,
     1695                   (long*) number_of_boundaries -> data,
     1696                   (double*) centroid_coordinates -> data,
     1697                   (double*) stage_centroid_values -> data,
     1698                   (double*) xmom_centroid_values -> data,
     1699                   (double*) ymom_centroid_values -> data,
     1700                   (double*) elevation_centroid_values -> data,
     1701                   (double*) vertex_coordinates -> data,
     1702                   (double*) stage_vertex_values -> data,
     1703                   (double*) xmom_vertex_values -> data,
     1704                   (double*) ymom_vertex_values -> data,
     1705                   (double*) elevation_vertex_values -> data,
     1706                   optimise_dry_cells);
    15801707  if (e == -1) {
    15811708    // Use error string set inside computational routine
    15821709    return NULL;
    1583   }                              
     1710  }                  
    15841711 
    15851712 
    … …  
    16091736  // Convert Python arguments to C
    16101737  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OO|i", argnames,
    1611                                    &Q, &Normal, &direction)) {
     1738                   &Q, &Normal, &direction)) {
    16121739    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: rotate could not parse input arguments");
    16131740    return NULL;
    … …  
    16541781// Computational function for flux computation
    16551782double _compute_fluxes_central(int number_of_elements,
    1656                                double timestep,
    1657                                double epsilon,
    1658                                double H0,
    1659                                double g,
    1660                                long* neighbours,
    1661                                long* neighbour_edges,
    1662                                double* normals,
    1663                                double* edgelengths,
    1664                                double* radii,
    1665                                double* areas,
    1666                                long* tri_full_flag,
    1667                                double* stage_edge_values,
    1668                                double* xmom_edge_values,
    1669                                double* ymom_edge_values,
    1670                                double* bed_edge_values,
    1671                                double* stage_boundary_values,
    1672                                double* xmom_boundary_values,
    1673                                double* ymom_boundary_values,
    1674                                double* stage_explicit_update,
    1675                                double* xmom_explicit_update,
    1676                                double* ymom_explicit_update,
    1677                                long* already_computed_flux,
    1678                                double* max_speed_array,
    1679                                int optimise_dry_cells) {
    1680                                
     1783                               double timestep,
     1784                               double epsilon,
     1785                               double H0,
     1786                               double g,
     1787                               long* neighbours,
     1788                               long* neighbour_edges,
     1789                               double* normals,
     1790                               double* edgelengths,
     1791                               double* radii,
     1792                               double* areas,
     1793                               long* tri_full_flag,
     1794                               double* stage_edge_values,
     1795                               double* xmom_edge_values,
     1796                               double* ymom_edge_values,
     1797                               double* bed_edge_values,
     1798                               double* stage_boundary_values,
     1799                               double* xmom_boundary_values,
     1800                               double* ymom_boundary_values,
     1801                               double* stage_explicit_update,
     1802                               double* xmom_explicit_update,
     1803                               double* ymom_explicit_update,
     1804                               long* already_computed_flux,
     1805                               double* max_speed_array,
     1806                               int optimise_dry_cells)
     1807{                   
    16811808  // Local variables
    1682   double max_speed, length, area, zl, zr;
     1809  double max_speed, length, inv_area, zl, zr;
    16831810  int k, i, m, n;
    16841811  int ki, nm=0, ki2; // Index shorthands
    … …  
    16871814  double ql[3], qr[3], edgeflux[3]; // Work array for summing up fluxes
    16881815
    1689   static long call=1; // Static local variable flagging already computed flux
    1690                                
    1691        
     1816  static long call = 1; // Static local variable flagging already computed flux
     1817                   
    16921818  // Start computation
    16931819  call++; // Flag 'id' of flux calculation for this timestep
    1694  
     1820
    16951821  // Set explicit_update to zero for all conserved_quantities.
    16961822  // This assumes compute_fluxes called before forcing terms
    1697   for (k=0; k<number_of_elements; k++) {
    1698     stage_explicit_update[k]=0.0;
    1699     xmom_explicit_update[k]=0.0;
    1700     ymom_explicit_update[k]=0.0; 
    1701   }
    1702 
     1823  memset((char*) stage_explicit_update, 0, number_of_elements*sizeof(double));
     1824  memset((char*) xmom_explicit_update, 0, number_of_elements*sizeof(double));
     1825  memset((char*) ymom_explicit_update, 0, number_of_elements*sizeof(double));   
     1826 
    17031827  // For all triangles
    1704   for (k=0; k<number_of_elements; k++) {
    1705    
     1828  for (k = 0; k < number_of_elements; k++)
     1829  { 
    17061830    // Loop through neighbours and compute edge flux for each 
    1707     for (i=0; i<3; i++) {
    1708       ki = k*3+i; // Linear index (triangle k, edge i)
     1831    for (i = 0; i < 3; i++)
     1832    {
     1833      ki = k*3 + i; // Linear index to edge i of triangle k
    17091834     
    17101835      if (already_computed_flux[ki] == call)
     1836      {
    17111837        // We've already computed the flux across this edge
    1712         continue;
    1713        
    1714        
     1838        continue;
     1839      }
     1840   
     1841      // Get left hand side values from triangle k, edge i
    17151842      ql[0] = stage_edge_values[ki];
    17161843      ql[1] = xmom_edge_values[ki];
    … …  
    17181845      zl = bed_edge_values[ki];
    17191846
    1720       // Quantities at neighbour on nearest face
     1847      // Get right hand side values either from neighbouring triangle
     1848      // or from boundary array (Quantities at neighbour on nearest face).
    17211849      n = neighbours[ki];
    1722       if (n < 0) {
     1850      if (n < 0)
     1851      {
    17231852        // Neighbour is a boundary condition
    1724         m = -n-1; // Convert negative flag to boundary index
    1725        
    1726         qr[0] = stage_boundary_values[m];
    1727         qr[1] = xmom_boundary_values[m];
    1728         qr[2] = ymom_boundary_values[m];
    1729         zr = zl; // Extend bed elevation to boundary
    1730       } else {
    1731         // Neighbour is a real element
    1732         m = neighbour_edges[ki];
    1733         nm = n*3+m; // Linear index (triangle n, edge m)
    1734        
    1735         qr[0] = stage_edge_values[nm];
    1736         qr[1] = xmom_edge_values[nm];
    1737         qr[2] = ymom_edge_values[nm];
    1738         zr = bed_edge_values[nm];
    1739       }
    1740      
    1741      
    1742       if (optimise_dry_cells) {     
    1743         // Check if flux calculation is necessary across this edge
    1744         // This check will exclude dry cells.
    1745         // This will also optimise cases where zl != zr as
    1746         // long as both are dry
    1747 
    1748         if ( fabs(ql[0] - zl) < epsilon &&
    1749              fabs(qr[0] - zr) < epsilon ) {
    1750           // Cell boundary is dry
    1751          
    1752           already_computed_flux[ki] = call; // #k Done 
    1753           if (n>=0)
    1754             already_computed_flux[nm] = call; // #n Done
    1755        
    1756           max_speed = 0.0;
    1757           continue;
    1758         }
    1759       }
    1760    
     1853        m = -n - 1; // Convert negative flag to boundary index
     1854   
     1855        qr[0] = stage_boundary_values[m];
     1856        qr[1] = xmom_boundary_values[m];
     1857        qr[2] = ymom_boundary_values[m];
     1858        zr = zl; // Extend bed elevation to boundary
     1859      }
     1860      else
     1861      {
     1862        // Neighbour is a real triangle
     1863        m = neighbour_edges[ki];
     1864        nm = n*3 + m; // Linear index (triangle n, edge m)
     1865       
     1866        qr[0] = stage_edge_values[nm];
     1867        qr[1] = xmom_edge_values[nm];
     1868        qr[2] = ymom_edge_values[nm];
     1869        zr = bed_edge_values[nm];
     1870      }
     1871     
     1872      // Now we have values for this edge - both from left and right side.
     1873     
     1874      if (optimise_dry_cells)
     1875      {     
     1876        // Check if flux calculation is necessary across this edge
     1877        // This check will exclude dry cells.
     1878        // This will also optimise cases where zl != zr as
     1879        // long as both are dry
     1880
     1881        if (fabs(ql[0] - zl) < epsilon &&
     1882            fabs(qr[0] - zr) < epsilon)
     1883        {
     1884          // Cell boundary is dry
     1885         
     1886          already_computed_flux[ki] = call; // #k Done 
     1887          if (n >= 0)
     1888          {
     1889            already_computed_flux[nm] = call; // #n Done
     1890          }
     1891       
     1892          max_speed = 0.0;
     1893          continue;
     1894        }
     1895      }
    17611896     
    17621897      // Outward pointing normal vector (domain.normals[k, 2*i:2*i+2])
    … …  
    17651900      // Edge flux computation (triangle k, edge i)
    17661901      _flux_function_central(ql, qr, zl, zr,
    1767                              normals[ki2], normals[ki2+1],
    1768                              epsilon, H0, g,
    1769                              edgeflux, &max_speed);
     1902                             normals[ki2], normals[ki2+1],
     1903                             epsilon, H0, g,
     1904                             edgeflux, &max_speed);
    17701905     
    17711906     
    … …  
    17861921     
    17871922      // Update neighbour n with same flux but reversed sign
    1788       if (n>=0) {
    1789         stage_explicit_update[n] += edgeflux[0];
    1790         xmom_explicit_update[n] += edgeflux[1];
    1791         ymom_explicit_update[n] += edgeflux[2];
    1792        
    1793         already_computed_flux[nm] = call; // #n Done
    1794       }
    1795 
     1923      if (n >= 0)
     1924      {
     1925        stage_explicit_update[n] += edgeflux[0];
     1926        xmom_explicit_update[n] += edgeflux[1];
     1927        ymom_explicit_update[n] += edgeflux[2];
     1928       
     1929        already_computed_flux[nm] = call; // #n Done
     1930      }
    17961931
    17971932      // Update timestep based on edge i and possibly neighbour n
    1798       if (tri_full_flag[k] == 1) {
    1799         if (max_speed > epsilon) {
    1800        
    1801           // Apply CFL condition for triangles joining this edge (triangle k and triangle n)
    1802          
    1803           // CFL for triangle k
    1804           timestep = min(timestep, radii[k]/max_speed);
    1805          
    1806           if (n>=0)
    1807             // Apply CFL condition for neigbour n (which is on the ith edge of triangle k)
    1808             timestep = min(timestep, radii[n]/max_speed);
    1809          
    1810           // Ted Rigby's suggested less conservative version
    1811           //if (n>=0) {             
    1812           //  timestep = min(timestep, (radii[k]+radii[n])/max_speed);
    1813           //} else {
    1814           //  timestep = min(timestep, radii[k]/max_speed);
    1815           // }
    1816         }
     1933      if (tri_full_flag[k] == 1)
     1934      {
     1935        if (max_speed > epsilon)
     1936        {
     1937          // Apply CFL condition for triangles joining this edge (triangle k and triangle n)
     1938         
     1939          // CFL for triangle k
     1940          timestep = min(timestep, radii[k]/max_speed);
     1941         
     1942          if (n >= 0)
     1943          {
     1944            // Apply CFL condition for neigbour n (which is on the ith edge of triangle k)
     1945            timestep = min(timestep, radii[n]/max_speed);
     1946          }
     1947
     1948          // Ted Rigby's suggested less conservative version
     1949          //if (n>=0) {         
     1950          //  timestep = min(timestep, (radii[k]+radii[n])/max_speed);
     1951          //} else {
     1952          //  timestep = min(timestep, radii[k]/max_speed);
     1953          // }
     1954        }
    18171955      }
    18181956     
    … …  
    18221960    // Normalise triangle k by area and store for when all conserved
    18231961    // quantities get updated
    1824     area = areas[k];
    1825     stage_explicit_update[k] /= area;
    1826     xmom_explicit_update[k] /= area;
    1827     ymom_explicit_update[k] /= area;
     1962    inv_area = 1.0/areas[k];
     1963    stage_explicit_update[k] *= inv_area;
     1964    xmom_explicit_update[k] *= inv_area;
     1965    ymom_explicit_update[k] *= inv_area;
    18281966   
    18291967   
    … …  
    18321970   
    18331971  } // End triangle k
    1834        
    1835                                
    1836                                
     1972             
    18371973  return timestep;
    18381974}
    … …  
    18692005
    18702006    PyObject
    1871         *domain,
    1872         *stage,
    1873         *xmom,
    1874         *ymom,
    1875         *bed;
     2007    *domain,
     2008    *stage,
     2009    *xmom,
     2010    *ymom,
     2011    *bed;
    18762012
    18772013    PyArrayObject
    1878         *neighbours,
    1879         *neighbour_edges,
    1880         *normals,
    1881         *edgelengths,
    1882         *radii,
    1883         *areas,
    1884         *tri_full_flag,
    1885         *stage_edge_values,
    1886         *xmom_edge_values,
    1887         *ymom_edge_values,
    1888         *bed_edge_values,
    1889         *stage_boundary_values,
    1890         *xmom_boundary_values,
    1891         *ymom_boundary_values,
    1892         *stage_explicit_update,
    1893         *xmom_explicit_update,
    1894         *ymom_explicit_update,
    1895         *already_computed_flux, //Tracks whether the flux across an edge has already been computed
    1896         *max_speed_array; //Keeps track of max speeds for each triangle
     2014    *neighbours,
     2015    *neighbour_edges,
     2016    *normals,
     2017    *edgelengths,
     2018    *radii,
     2019    *areas,
     2020    *tri_full_flag,
     2021    *stage_edge_values,
     2022    *xmom_edge_values,
     2023    *ymom_edge_values,
     2024    *bed_edge_values,
     2025    *stage_boundary_values,
     2026    *xmom_boundary_values,
     2027    *ymom_boundary_values,
     2028    *stage_explicit_update,
     2029    *xmom_explicit_update,
     2030    *ymom_explicit_update,
     2031    *already_computed_flux, //Tracks whether the flux across an edge has already been computed
     2032    *max_speed_array; //Keeps track of max speeds for each triangle
    18972033
    18982034   
    … …  
    19022038    // Convert Python arguments to C
    19032039    if (!PyArg_ParseTuple(args, "dOOOO", &timestep, &domain, &stage, &xmom, &ymom, &bed )) {
    1904         PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
    1905         return NULL;
     2040    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
     2041    return NULL;
    19062042    }
    19072043
    … …  
    19402076  // the explicit update arrays
    19412077  timestep = _compute_fluxes_central(number_of_elements,
    1942                                      timestep,
    1943                                      epsilon,
    1944                                      H0,
    1945                                      g,
    1946                                      (long*) neighbours -> data,
    1947                                      (long*) neighbour_edges -> data,
    1948                                      (double*) normals -> data,
    1949                                      (double*) edgelengths -> data,
    1950                                      (double*) radii -> data,
    1951                                      (double*) areas -> data,
    1952                                      (long*) tri_full_flag -> data,
    1953                                      (double*) stage_edge_values -> data,
    1954                                      (double*) xmom_edge_values -> data,
    1955                                      (double*) ymom_edge_values -> data,
    1956                                      (double*) bed_edge_values -> data,
    1957                                      (double*) stage_boundary_values -> data,
    1958                                      (double*) xmom_boundary_values -> data,
    1959                                      (double*) ymom_boundary_values -> data,
    1960                                      (double*) stage_explicit_update -> data,
    1961                                      (double*) xmom_explicit_update -> data,
    1962                                      (double*) ymom_explicit_update -> data,
    1963                                      (long*) already_computed_flux -> data,
    1964                                      (double*) max_speed_array -> data,
    1965                                      optimise_dry_cells);
     2078                     timestep,
     2079                     epsilon,
     2080                     H0,
     2081                     g,
     2082                     (long*) neighbours -> data,
     2083                     (long*) neighbour_edges -> data,
     2084                     (double*) normals -> data,
     2085                     (double*) edgelengths -> data,
     2086                     (double*) radii -> data,
     2087                     (double*) areas -> data,
     2088                     (long*) tri_full_flag -> data,
     2089                     (double*) stage_edge_values -> data,
     2090                     (double*) xmom_edge_values -> data,
     2091                     (double*) ymom_edge_values -> data,
     2092                     (double*) bed_edge_values -> data,
     2093                     (double*) stage_boundary_values -> data,
     2094                     (double*) xmom_boundary_values -> data,
     2095                     (double*) ymom_boundary_values -> data,
     2096                     (double*) stage_explicit_update -> data,
     2097                     (double*) xmom_explicit_update -> data,
     2098                     (double*) ymom_explicit_update -> data,
     2099                     (long*) already_computed_flux -> data,
     2100                     (double*) max_speed_array -> data,
     2101                     optimise_dry_cells);
    19662102
    19672103  Py_DECREF(neighbours);
    … …  
    20132149    domain.timestep = compute_fluxes(timestep,
    20142150                                     domain.epsilon,
    2015                                      domain.H0,
     2151                     domain.H0,
    20162152                                     domain.g,
    20172153                                     domain.neighbours,
    … …  
    20332169                                     Ymom.explicit_update,
    20342170                                     already_computed_flux,
    2035                                      optimise_dry_cells)                                     
     2171                     optimise_dry_cells)                     
    20362172
    20372173
    … …  
    20662202  // Convert Python arguments to C
    20672203  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOi",
    2068                         &timestep,
    2069                         &epsilon,
    2070                         &H0,
    2071                         &g,
    2072                         &neighbours,
    2073                         &neighbour_edges,
    2074                         &normals,
    2075                         &edgelengths, &radii, &areas,
    2076                         &tri_full_flag,
    2077                         &stage_edge_values,
    2078                         &xmom_edge_values,
    2079                         &ymom_edge_values,
    2080                         &bed_edge_values,
    2081                         &stage_boundary_values,
    2082                         &xmom_boundary_values,
    2083                         &ymom_boundary_values,
    2084                         &stage_explicit_update,
    2085                         &xmom_explicit_update,
    2086                         &ymom_explicit_update,
    2087                         &already_computed_flux,
    2088                         &max_speed_array,
    2089                         &optimise_dry_cells)) {
     2204            &timestep,
     2205            &epsilon,
     2206            &H0,
     2207            &g,
     2208            &neighbours,
     2209            &neighbour_edges,
     2210            &normals,
     2211            &edgelengths, &radii, &areas,
     2212            &tri_full_flag,
     2213            &stage_edge_values,
     2214            &xmom_edge_values,
     2215            &ymom_edge_values,
     2216            &bed_edge_values,
     2217            &stage_boundary_values,
     2218            &xmom_boundary_values,
     2219            &ymom_boundary_values,
     2220            &stage_explicit_update,
     2221            &xmom_explicit_update,
     2222            &ymom_explicit_update,
     2223            &already_computed_flux,
     2224            &max_speed_array,
     2225            &optimise_dry_cells)) {
    20902226    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
    20912227    return NULL;
    … …  
    21182254  // the explicit update arrays
    21192255  timestep = _compute_fluxes_central(number_of_elements,
    2120                                      timestep,
    2121                                      epsilon,
    2122                                      H0,
    2123                                      g,
    2124                                      (long*) neighbours -> data,
    2125                                      (long*) neighbour_edges -> data,
    2126                                      (double*) normals -> data,
    2127                                      (double*) edgelengths -> data,
    2128                                      (double*) radii -> data,
    2129                                      (double*) areas -> data,
    2130                                      (long*) tri_full_flag -> data,
    2131                                      (double*) stage_edge_values -> data,
    2132                                      (double*) xmom_edge_values -> data,
    2133                                      (double*) ymom_edge_values -> data,
    2134                                      (double*) bed_edge_values -> data,
    2135                                      (double*) stage_boundary_values -> data,
    2136                                      (double*) xmom_boundary_values -> data,
    2137                                      (double*) ymom_boundary_values -> data,
    2138                                      (double*) stage_explicit_update -> data,
    2139                                      (double*) xmom_explicit_update -> data,
    2140                                      (double*) ymom_explicit_update -> data,
    2141                                      (long*) already_computed_flux -> data,
    2142                                      (double*) max_speed_array -> data,
    2143                                      optimise_dry_cells);
     2256                     timestep,
     2257                     epsilon,
     2258                     H0,
     2259                     g,
     2260                     (long*) neighbours -> data,
     2261                     (long*) neighbour_edges -> data,
     2262                     (double*) normals -> data,
     2263                     (double*) edgelengths -> data,
     2264                     (double*) radii -> data,
     2265                     (double*) areas -> data,
     2266                     (long*) tri_full_flag -> data,
     2267                     (double*) stage_edge_values -> data,
     2268                     (double*) xmom_edge_values -> data,
     2269                     (double*) ymom_edge_values -> data,
     2270                     (double*) bed_edge_values -> data,
     2271                     (double*) stage_boundary_values -> data,
     2272                     (double*) xmom_boundary_values -> data,
     2273                     (double*) ymom_boundary_values -> data,
     2274                     (double*) stage_explicit_update -> data,
     2275                     (double*) xmom_explicit_update -> data,
     2276                     (double*) ymom_explicit_update -> data,
     2277                     (long*) already_computed_flux -> data,
     2278                     (double*) max_speed_array -> data,
     2279                     optimise_dry_cells);
    21442280 
    21452281  // Return updated flux timestep
    … …  
    22282364  // Convert Python arguments to C
    22292365  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOO",
    2230                         &timestep,
    2231                         &epsilon,
    2232                         &H0,
    2233                         &g,
    2234                         &neighbours,
    2235                         &neighbour_edges,
    2236                         &normals,
    2237                         &edgelengths, &radii, &areas,
    2238                         &tri_full_flag,
    2239                         &stage_edge_values,
    2240                         &xmom_edge_values,
    2241                         &ymom_edge_values,
    2242                         &bed_edge_values,
    2243                         &stage_boundary_values,
    2244                         &xmom_boundary_values,
    2245                         &ymom_boundary_values,
    2246                         &stage_explicit_update,
    2247                         &xmom_explicit_update,
    2248                         &ymom_explicit_update,
    2249                         &already_computed_flux)) {
     2366            &timestep,
     2367            &epsilon,
     2368            &H0,
     2369            &g,
     2370            &neighbours,
     2371            &neighbour_edges,
     2372            &normals,
     2373            &edgelengths, &radii, &areas,
     2374            &tri_full_flag,
     2375            &stage_edge_values,
     2376            &xmom_edge_values,
     2377            &ymom_edge_values,
     2378            &bed_edge_values,
     2379            &stage_boundary_values,
     2380            &xmom_boundary_values,
     2381            &ymom_boundary_values,
     2382            &stage_explicit_update,
     2383            &xmom_explicit_update,
     2384            &ymom_explicit_update,
     2385            &already_computed_flux)) {
    22502386    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
    22512387    return NULL;
    … …  
    22652401      ki = k*3+i;
    22662402      if (((long *) already_computed_flux->data)[ki]==call)//we've already computed the flux across this edge
    2267         continue;
     2403    continue;
    22682404      ql[0] = ((double *) stage_edge_values -> data)[ki];
    22692405      ql[1] = ((double *) xmom_edge_values -> data)[ki];
    … …  
    22742410      n = ((long *) neighbours -> data)[ki];
    22752411      if (n < 0) {
    2276         m = -n-1; //Convert negative flag to index
    2277         qr[0] = ((double *) stage_boundary_values -> data)[m];
    2278         qr[1] = ((double *) xmom_boundary_values -> data)[m];
    2279         qr[2] = ((double *) ymom_boundary_values -> data)[m];
    2280         zr = zl; //Extend bed elevation to boundary
     2412    m = -n-1; //Convert negative flag to index
     2413    qr[0] = ((double *) stage_boundary_values -> data)[m];
     2414    qr[1] = ((double *) xmom_boundary_values -> data)[m];
     2415    qr[2] = ((double *) ymom_boundary_values -> data)[m];
     2416    zr = zl; //Extend bed elevation to boundary
    22812417      } else {
    2282         m = ((long *) neighbour_edges -> data)[ki];
    2283         nm = n*3+m;
    2284         qr[0] = ((double *) stage_edge_values -> data)[nm];
    2285         qr[1] = ((double *) xmom_edge_values -> data)[nm];
    2286         qr[2] = ((double *) ymom_edge_values -> data)[nm];
    2287         zr =    ((double *) bed_edge_values -> data)[nm];
     2418    m = ((long *) neighbour_edges -> data)[ki];
     2419    nm = n*3+m;
     2420    qr[0] = ((double *) stage_edge_values -> data)[nm];
     2421    qr[1] = ((double *) xmom_edge_values -> data)[nm];
     2422    qr[2] = ((double *) ymom_edge_values -> data)[nm];
     2423    zr =    ((double *) bed_edge_values -> data)[nm];
    22882424      }
    22892425      // Outward pointing normal vector
    … …  
    22942430      //Edge flux computation
    22952431      flux_function_kinetic(ql, qr, zl, zr,
    2296                     normal[0], normal[1],
    2297                     epsilon, H0, g,
    2298                     edgeflux, &max_speed);
     2432            normal[0], normal[1],
     2433            epsilon, H0, g,
     2434            edgeflux, &max_speed);
    22992435      //update triangle k
    23002436      ((long *) already_computed_flux->data)[ki]=call;
    … …  
    23042440      //update the neighbour n
    23052441      if (n>=0){
    2306         ((long *) already_computed_flux->data)[nm]=call;
    2307         ((double *) stage_explicit_update -> data)[n] += edgeflux[0]*((double *) edgelengths -> data)[nm];
    2308         ((double *) xmom_explicit_update -> data)[n] += edgeflux[1]*((double *) edgelengths -> data)[nm];
    2309         ((double *) ymom_explicit_update -> data)[n] += edgeflux[2]*((double *) edgelengths -> data)[nm];
     2442    ((long *) already_computed_flux->data)[nm]=call;
     2443    ((double *) stage_explicit_update -> data)[n] += edgeflux[0]*((double *) edgelengths -> data)[nm];
     2444    ((double *) xmom_explicit_update -> data)[n] += edgeflux[1]*((double *) edgelengths -> data)[nm];
     2445    ((double *) ymom_explicit_update -> data)[n] += edgeflux[2]*((double *) edgelengths -> data)[nm];
    23102446      }
    23112447      ///for (j=0; j<3; j++) {
    2312         ///flux[j] -= edgeflux[j]*((double *) edgelengths -> data)[ki];
    2313         ///}
    2314         //Update timestep
    2315         //timestep = min(timestep, domain.radii[k]/max_speed)
    2316         //FIXME: SR Add parameter for CFL condition
     2448    ///flux[j] -= edgeflux[j]*((double *) edgelengths -> data)[ki];
     2449    ///}
     2450    //Update timestep
     2451    //timestep = min(timestep, domain.radii[k]/max_speed)
     2452    //FIXME: SR Add parameter for CFL condition
    23172453    if ( ((long *) tri_full_flag -> data)[k] == 1) {
    2318             if (max_speed > epsilon) {
    2319                 timestep = min(timestep, ((double *) radii -> data)[k]/max_speed);
    2320                 //maxspeed in flux_function is calculated as max(|u+a|,|u-a|)
    2321                 if (n>=0)
    2322                     timestep = min(timestep, ((double *) radii -> data)[n]/max_speed);
    2323             }
     2454        if (max_speed > epsilon) {
     2455            timestep = min(timestep, ((double *) radii -> data)[k]/max_speed);
     2456            //maxspeed in flux_function is calculated as max(|u+a|,|u-a|)
     2457            if (n>=0)
     2458                timestep = min(timestep, ((double *) radii -> data)[n]/max_speed);
     2459        }
    23242460    }
    23252461    } // end for i
    … …  
    23502486  // Convert Python arguments to C
    23512487  if (!PyArg_ParseTuple(args, "dddOOOO",
    2352                         &minimum_allowed_height,
    2353                         &maximum_allowed_speed,
    2354                         &epsilon,
    2355                         &wc, &zc, &xmomc, &ymomc)) {
     2488            &minimum_allowed_height,
     2489            &maximum_allowed_speed,
     2490            &epsilon,
     2491            &wc, &zc, &xmomc, &ymomc)) {
    23562492    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: protect could not parse input arguments");
    23572493    return NULL;
    … …  
    23612497
    23622498  _protect(N,
    2363            minimum_allowed_height,
    2364            maximum_allowed_speed,
    2365            epsilon,
    2366            (double*) wc -> data,
    2367            (double*) zc -> data,
    2368            (double*) xmomc -> data,
    2369            (double*) ymomc -> data);
     2499       minimum_allowed_height,
     2500       maximum_allowed_speed,
     2501       epsilon,
     2502       (double*) wc -> data,
     2503       (double*) zc -> data,
     2504       (double*) xmomc -> data,
     2505       (double*) ymomc -> data);
    23702506
    23712507  return Py_BuildValue("");
    … …  
    24092545  // Convert Python arguments to C
    24102546  if (!PyArg_ParseTuple(args, "OOOOOOOOOO",
    2411                         &domain,
    2412                         &wc, &zc,
    2413                         &wv, &zv, &hvbar,
    2414                         &xmomc, &ymomc, &xmomv, &ymomv)) {
     2547            &domain,
     2548            &wc, &zc,
     2549            &wv, &zv, &hvbar,
     2550            &xmomc, &ymomc, &xmomv, &ymomv)) {
    24152551    PyErr_SetString(PyExc_RuntimeError,
    2416                     "shallow_water_ext.c: balance_deep_and_shallow could not parse input arguments");
     2552            "shallow_water_ext.c: balance_deep_and_shallow could not parse input arguments");
    24172553    return NULL;
    24182554  } 
    2419          
    2420          
     2555     
     2556     
    24212557  // FIXME (Ole): I tested this without GetAttrString and got time down
    24222558  // marginally from 4.0s to 3.8s. Consider passing everything in
    … …  
    24632599  N = wc -> dimensions[0];
    24642600  _balance_deep_and_shallow(N,
    2465                             (double*) wc -> data,
    2466                             (double*) zc -> data,
    2467                             (double*) wv -> data,
    2468                             (double*) zv -> data,
    2469                             (double*) hvbar -> data,
    2470                             (double*) xmomc -> data,
    2471                             (double*) ymomc -> data,
    2472                             (double*) xmomv -> data,
    2473                             (double*) ymomv -> data,
    2474                             H0,
    2475                             (int) tight_slope_limiters,
    2476                             (int) use_centroid_velocities,                         
    2477                             alpha_balance);
     2601                (double*) wc -> data,
     2602                (double*) zc -> data,
     2603                (double*) wv -> data,
     2604                (double*) zv -> data,
     2605                (double*) hvbar -> data,
     2606                (double*) xmomc -> data,
     2607                (double*) ymomc -> data,
     2608                (double*) xmomv -> data,
     2609                (double*) ymomv -> data,
     2610                H0,
     2611                (int) tight_slope_limiters,
     2612                (int) use_centroid_velocities,             
     2613                alpha_balance);
    24782614
    24792615
    … …  
    25032639  // Convert Python arguments to C
    25042640  if (!PyArg_ParseTuple(args, "OOOOd",
    2505                         &xmom_update,
    2506                         &ymom_update,
    2507                         &s_vec, &phi_vec,
    2508                         &cw)) {
     2641            &xmom_update,
     2642            &ymom_update,
     2643            &s_vec, &phi_vec,
     2644            &cw)) {
    25092645    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: assign_windfield_values could not parse input arguments");
    25102646    return NULL;
    25112647  } 
    2512                        
     2648           
    25132649
    25142650  N = xmom_update -> dimensions[0];
    25152651
    25162652  _assign_wind_field_values(N,
    2517            (double*) xmom_update -> data,
    2518            (double*) ymom_update -> data,
    2519            (double*) s_vec -> data,
    2520            (double*) phi_vec -> data,
    2521            cw);
     2653       (double*) xmom_update -> data,
     2654       (double*) ymom_update -> data,
     2655       (double*) s_vec -> data,
     2656       (double*) phi_vec -> data,
     2657       cw);
    25222658
    25232659  return Py_BuildValue("");
Note: See TracChangeset for help on using the changeset viewer.