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

Back-merge from Numeric trunk to numpy branch.

Location:
branches/numpy/anuga/shallow_water
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • branches/numpy/anuga/shallow_water/data_manager.py

    r6689 r6902  
    35753575    other_quantities.remove('z')
    35763576    other_quantities.remove('volumes')
    3577     #try:
    3578     #    other_quantities.remove('stage_range')
    3579     #    other_quantities.remove('xmomentum_range')
    3580     #    other_quantities.remove('ymomentum_range')
    3581     #    other_quantities.remove('elevation_range')
    3582     #except:
    3583     #    pass
     3577    try:
     3578        other_quantities.remove('stage_range')
     3579        other_quantities.remove('xmomentum_range')
     3580        other_quantities.remove('ymomentum_range')
     3581        other_quantities.remove('elevation_range')
     3582    except:
     3583        pass
    35843584
    35853585    conserved_quantities.remove('time')
     
    53385338        j += 1
    53395339
    5340     #if verbose: sww.verbose_quantities(outfile)
     5340    if verbose: sww.verbose_quantities(outfile)
    53415341
    53425342    outfile.close()
     
    58555855    # @param smoothing True if smoothing is to be used.
    58565856    # @param order
    5857     # @param sww_precision Data type of the quantitiy written (netcdf constant)
     5857    # @param sww_precision Data type of the quantity written (netcdf constant)
    58585858    # @param verbose True if this function is to be verbose.
    58595859    # @note If 'times' is a list, the info will be made relative.
     
    59365936                               ('number_of_points',))
    59375937        q = 'elevation'
    5938         #outfile.createVariable(q + Write_sww.RANGE, sww_precision,
    5939         #                       ('numbers_in_range',))
    5940 
    5941 
    5942         ## Initialise ranges with small and large sentinels.
    5943         ## If this was in pure Python we could have used None sensibly
    5944         #outfile.variables[q+Write_sww.RANGE][0] = max_float  # Min
    5945         #outfile.variables[q+Write_sww.RANGE][1] = -max_float # Max
     5938        outfile.createVariable(q + Write_sww.RANGE, sww_precision,
     5939                               ('numbers_in_range',))
     5940
     5941        # Initialise ranges with small and large sentinels.
     5942        # If this was in pure Python we could have used None sensibly
     5943        outfile.variables[q+Write_sww.RANGE][0] = max_float  # Min
     5944        outfile.variables[q+Write_sww.RANGE][1] = -max_float # Max
    59465945
    59475946        # FIXME: Backwards compatibility
     
    59585957            outfile.createVariable(q, sww_precision, ('number_of_timesteps',
    59595958                                                      'number_of_points'))
    5960             #outfile.createVariable(q + Write_sww.RANGE, sww_precision,
    5961             #                       ('numbers_in_range',))
     5959            outfile.createVariable(q + Write_sww.RANGE, sww_precision,
     5960                                   ('numbers_in_range',))
    59625961
    59635962            # Initialise ranges with small and large sentinels.
    59645963            # If this was in pure Python we could have used None sensibly
    5965             #outfile.variables[q+Write_sww.RANGE][0] = max_float  # Min
    5966             #outfile.variables[q+Write_sww.RANGE][1] = -max_float # Max
     5964            outfile.variables[q+Write_sww.RANGE][0] = max_float  # Min
     5965            outfile.variables[q+Write_sww.RANGE][1] = -max_float # Max
    59675966
    59685967        if isinstance(times, (list, num.ndarray)):
     
    60736072        q = 'elevation'
    60746073        # This updates the _range values
    6075         #outfile.variables[q + Write_sww.RANGE][0] = min(elevation)
    6076         #outfile.variables[q + Write_sww.RANGE][1] = max(elevation)
     6074        outfile.variables[q + Write_sww.RANGE][0] = num.min(elevation)
     6075        outfile.variables[q + Write_sww.RANGE][1] = num.max(elevation)
    60776076
    60786077
     
    61266125                                q_values.astype(sww_precision)
    61276126       
    6128         #        # This updates the _range values
    6129         #        q_range = outfile.variables[q + Write_sww.RANGE][:]
    6130         #        q_values_min = min(q_values)
    6131         #        if q_values_min < q_range[0]:
    6132         #            outfile.variables[q + Write_sww.RANGE][0] = q_values_min
    6133         #        q_values_max = max(q_values)
    6134         #        if q_values_max > q_range[1]:
    6135         #            outfile.variables[q + Write_sww.RANGE][1] = q_values_max
     6127                # This updates the _range values
     6128                q_range = outfile.variables[q + Write_sww.RANGE][:]
     6129                q_values_min = num.min(q_values)
     6130                if q_values_min < q_range[0]:
     6131                    outfile.variables[q + Write_sww.RANGE][0] = q_values_min
     6132                q_values_max = num.max(q_values)
     6133                if q_values_max > q_range[1]:
     6134                    outfile.variables[q + Write_sww.RANGE][1] = q_values_max
    61366135
    61376136    ##
    61386137    # @brief Print the quantities in the underlying file.
    61396138    # @param outfile UNUSED.
    6140     #def verbose_quantities(self, outfile):
    6141     #    print '------------------------------------------------'
    6142     #    print 'More Statistics:'
    6143     #    for q in Write_sww.sww_quantities:
    6144     #        print '  %s in [%f, %f]' % (q,
    6145     #                                    outfile.variables[q+Write_sww.RANGE][0],
    6146     #                                    outfile.variables[q+Write_sww.RANGE][1])
    6147     #    print '------------------------------------------------'
     6139    def verbose_quantities(self, outfile):
     6140        print '------------------------------------------------'
     6141        print 'More Statistics:'
     6142        for q in Write_sww.sww_quantities:
     6143            print '  %s in [%f, %f]' % (q,
     6144                                        outfile.variables[q+Write_sww.RANGE][0],
     6145                                        outfile.variables[q+Write_sww.RANGE][1])
     6146        print '------------------------------------------------'
    61486147
    61496148
     
    64966495
    64976496                # This updates the _range values
    6498                 #q_range = outfile.variables[q + Write_sts.RANGE][:]
    6499                 #q_values_min = min(q_values)
    6500                 #if q_values_min < q_range[0]:
    6501                 #    outfile.variables[q + Write_sts.RANGE][0] = q_values_min
    6502                 #q_values_max = max(q_values)
    6503                 #if q_values_max > q_range[1]:
    6504                 #    outfile.variables[q + Write_sts.RANGE][1] = q_values_max
     6497                q_range = outfile.variables[q + Write_sts.RANGE][:]
     6498                q_values_min = num.min(q_values)
     6499                if q_values_min < q_range[0]:
     6500                    outfile.variables[q + Write_sts.RANGE][0] = q_values_min
     6501                q_values_max = num.max(q_values)
     6502                if q_values_max > q_range[1]:
     6503                    outfile.variables[q + Write_sts.RANGE][1] = q_values_max
    65056504
    65066505
     
    66736672
    66746673
    6675 # FIXME (DSG): Add unit test, make general, not just 2 files,
    6676 # but any number of files.
    66776674##
    66786675# @brief Copy a file to a directory, and optionally append another file to it.
     
    66806677# @param filename Path to file to copy to directory 'dir_name'.
    66816678# @param filename2 Optional path to file to append to copied file.
    6682 def copy_code_files(dir_name, filename1, filename2=None):
     6679# @param verbose True if this function is to be verbose.
     6680# @note Allow filenames to be either a string or sequence of strings.
     6681def copy_code_files(dir_name, filename1, filename2=None, verbose=False):
    66836682    """Copies "filename1" and "filename2" to "dir_name".
    6684     Very useful for information management
    6685     filename1 and filename2 are both absolute pathnames
    6686     """
    6687 
     6683
     6684    Each 'filename' may be a string or list of filename strings.
     6685
     6686    Filenames must be absolute pathnames
     6687    """
     6688
     6689    ##
     6690    # @brief copies a file or sequence to destination directory.
     6691    # @param dest The destination directory to copy to.
     6692    # @param file A filename string or sequence of filename strings.
     6693    def copy_file_or_sequence(dest, file):
     6694        if hasattr(file, '__iter__'):
     6695            for f in file:
     6696                shutil.copy(f, dir_name)
     6697                if verbose:
     6698                    print 'File %s copied' % (f)
     6699        else:
     6700            shutil.copy(file, dir_name)
     6701            if verbose:
     6702                print 'File %s copied' % (file)
     6703
     6704    # check we have a destination directory, create if necessary
    66886705    if access(dir_name, F_OK) == 0:
    6689         print 'Make directory %s' % dir_name
    6690         mkdir (dir_name,0777)
    6691 
    6692     shutil.copy(filename1, dir_name + sep + basename(filename1))
    6693 
    6694     if filename2 != None:
    6695         shutil.copy(filename2, dir_name + sep + basename(filename2))
    6696         print 'Files %s and %s copied' %(filename1, filename2)
    6697     else:
    6698         print 'File %s copied' %(filename1)
     6706        if verbose:
     6707            print 'Make directory %s' % dir_name
     6708        mkdir(dir_name, 0777)
     6709
     6710    copy_file_or_sequence(dir_name, filename1)
     6711
     6712    if not filename2 is None:
     6713        copy_file_or_sequence(dir_name, filename2)
    66996714
    67006715
  • 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("");
  • branches/numpy/anuga/shallow_water/test_data_manager.py

    r6689 r6902  
    1515from struct import pack, unpack
    1616from sets import ImmutableSet
     17import shutil
    1718
    1819from anuga.shallow_water import *
     
    4748       
    4849        self.verbose = Test_Data_Manager.verbose
    49         #Create basic mesh
     50        # Create basic mesh
    5051        points, vertices, boundary = rectangular(2, 2)
    5152
    52         #Create shallow water domain
     53        # Create shallow water domain
    5354        domain = Domain(points, vertices, boundary)
    5455        domain.default_order = 2
    5556
    56         #Set some field values
     57        # Set some field values
    5758        domain.set_quantity('elevation', lambda x,y: -x)
    5859        domain.set_quantity('friction', 0.03)
     
    252253
    253254        # Get the variables
    254         #range = fid.variables['stage_range'][:]
     255        range = fid.variables['stage_range'][:]
    255256        #print range
    256         #assert num.allclose(range,[-0.93519, 0.15]) or\
    257         #       num.allclose(range,[-0.9352743, 0.15]) or\
    258         #       num.allclose(range,[-0.93522203, 0.15000001]) # Old slope limiters
    259        
    260         #range = fid.variables['xmomentum_range'][:]
     257        assert num.allclose(range,[-0.93519, 0.15]) or\
     258               num.allclose(range,[-0.9352743, 0.15]) or\
     259               num.allclose(range,[-0.93522203, 0.15000001]) # Old slope limiters
     260       
     261        range = fid.variables['xmomentum_range'][:]
    261262        ##print range
    262         #assert num.allclose(range,[0,0.4695096]) or \
    263         #       num.allclose(range,[0,0.47790655]) or\
    264         #       num.allclose(range,[0,0.46941957]) or\
    265         #       num.allclose(range,[0,0.47769409])
    266 
    267        
    268         #range = fid.variables['ymomentum_range'][:]
     263        assert num.allclose(range,[0,0.4695096]) or \
     264               num.allclose(range,[0,0.47790655]) or\
     265               num.allclose(range,[0,0.46941957]) or\
     266               num.allclose(range,[0,0.47769409])
     267
     268       
     269        range = fid.variables['ymomentum_range'][:]
    269270        ##print range
    270         #assert num.allclose(range,[0,0.02174380]) or\
    271         #       num.allclose(range,[0,0.02174439]) or\
    272         #       num.allclose(range,[0,0.02283983]) or\
    273         #       num.allclose(range,[0,0.0217342]) or\
    274         #       num.allclose(range,[0,0.0227564]) # Old slope limiters
     271        assert num.allclose(range,[0,0.02174380]) or\
     272               num.allclose(range,[0,0.02174439]) or\
     273               num.allclose(range,[0,0.02283983]) or\
     274               num.allclose(range,[0,0.0217342]) or\
     275               num.allclose(range,[0,0.0227564]) # Old slope limiters
    275276       
    276277        fid.close()
     
    12321233        from Scientific.IO.NetCDF import NetCDFFile
    12331234
    1234         #Setup
     1235        # Setup
    12351236        self.domain.set_name('datatest')
    12361237
     
    12451246        self.domain.set_quantity('stage', 1.0)
    12461247
    1247         self.domain.geo_reference = Geo_reference(56,308500,6189000)
     1248        self.domain.geo_reference = Geo_reference(56, 308500, 6189000)
    12481249
    12491250        sww = get_dataobject(self.domain)
     
    55635564            quantities_init[0].append(this_ha) # HA
    55645565            quantities_init[1].append(this_ua) # UA
    5565             quantities_init[2].append(this_va) # VA
     5566            quantities_init[2].append(this_va) # VA 
    55665567               
    55675568        file_handle, base_name = tempfile.mkstemp("")
     
    61046105                        #print 'writing', time, point_i, q_time[time, point_i]
    61056106                        f.write(pack('f', q_time[time, point_i]))
    6106 
    61076107            f.close()
    61086108
     
    62156215
    62166216        msg='Incorrect gauge depths returned'
    6217         assert num.allclose(elevation,-depth),msg
     6217        assert num.allclose(elevation, -depth),msg
    62186218        msg='incorrect gauge height time series returned'
    6219         assert num.allclose(stage,ha)
     6219        assert num.allclose(stage, ha)
    62206220        msg='incorrect gauge ua time series returned'
    6221         assert num.allclose(xvelocity,ua)
     6221        assert num.allclose(xvelocity, ua)
    62226222        msg='incorrect gauge va time series returned'
    62236223        assert num.allclose(yvelocity, -va) # South is positive in MUX
     
    63976397                if j == 2: assert num.allclose(data[i][:parameters_index], -va0[permutation[i], :])
    63986398       
     6399        self.delete_mux(filesI)
    63996400
    64006401
     
    64076408        wrong values Win32
    64086409
    6409         This test does not pass on Windows but test_read_mux_platform_problem1 does
     6410        This test does not pass on Windows but test_read_mux_platform_problem1
     6411        does
    64106412        """
    64116413       
     
    64226424        times_ref = num.arange(0, time_step_count*time_step, time_step)
    64236425       
    6424         lat_long_points = [(-21.5,114.5), (-21,114.5), (-21.5,115), (-21.,115.), (-22., 117.)]
     6426        lat_long_points = [(-21.5,114.5), (-21,114.5), (-21.5,115),
     6427                           (-21.,115.), (-22., 117.)]
    64256428        n = len(lat_long_points)
    64266429       
     
    64666469             for j in range(0, first_tstep[i]-1) + range(last_tstep[i], time_step_count):
    64676470                 # For timesteps before and after recording range
    6468                  ha1[i][j] = ua1[i][j] = va1[i][j] = 0.0                                  
     6471                 ha1[i][j] = ua1[i][j] = va1[i][j] = 0.0
    64696472
    64706473
     
    64746477        #print 'va', va1[0,:]
    64756478       
    6476         # Write second mux file to be combined by urs2sts                                            
     6479        # Write second mux file to be combined by urs2sts
    64776480        base_nameII, filesII = self.write_mux2(lat_long_points,
    64786481                                               time_step_count, time_step,
     
    66056608
    66066609            f.close()
    6607 
    6608 
    6609 
    6610 
    6611 
    6612 
    6613        
    66146610                                               
    66156611        # Create ordering file
    66166612        permutation = ensure_numeric([4,0,2])
    66176613
    6618         _, ordering_filename = tempfile.mkstemp('')
    6619         order_fid = open(ordering_filename, 'w') 
    6620         order_fid.write('index, longitude, latitude\n')
    6621         for index in permutation:
    6622             order_fid.write('%d, %f, %f\n' %(index,
    6623                                              lat_long_points[index][1],
    6624                                              lat_long_points[index][0]))
    6625         order_fid.close()
    6626        
    6627        
    6628 
     6614       #  _, ordering_filename = tempfile.mkstemp('')
     6615#         order_fid = open(ordering_filename, 'w') 
     6616#         order_fid.write('index, longitude, latitude\n')
     6617#         for index in permutation:
     6618#             order_fid.write('%d, %f, %f\n' %(index,
     6619#                                              lat_long_points[index][1],
     6620#                                              lat_long_points[index][0]))
     6621#         order_fid.close()
     6622       
    66296623        # -------------------------------------
    66306624        # Now read files back and check values
     
    66386632        for j, file in enumerate(filesII):
    66396633            # Read stage, u, v enumerated as j
    6640 
    6641            
    66426634            #print 'Reading', j, file
    66436635            data = read_mux2(1, [file], weights, file_params, permutation, verbose)
     
    66456637            #print 'Data received by Python'
    66466638            #print data[1][8]
    6647 
    6648            
    66496639            number_of_selected_stations = data.shape[0]
    66506640
     
    66596649                #print i, parameters_index
    66606650                #print quantity[i][:]
    6661 
    6662                
    66636651                if j == 0: assert num.allclose(data[i][:parameters_index], ha1[permutation[i], :])
    66646652                if j == 1: assert num.allclose(data[i][:parameters_index], ua1[permutation[i], :])
     
    66696657                    #print j, i
    66706658                    #print 'Input'
    6671                     #print 'u', ua1[permutation[i], 8]                                        
     6659                    #print 'u', ua1[permutation[i], 8]       
    66726660                    #print 'v', va1[permutation[i], 8]
    66736661               
    66746662                    #print 'Output'
    6675                     #print 'v ', data[i][:parameters_index][8]                   
     6663                    #print 'v ', data[i][:parameters_index][8] 
     6664
     6665                    # South is positive in MUX
     6666                    #print "data[i][:parameters_index]", data[i][:parameters_index]
     6667                    #print "-va1[permutation[i], :]", -va1[permutation[i], :]
     6668                    assert num.allclose(data[i][:parameters_index], -va1[permutation[i], :])
     6669       
     6670        self.delete_mux(filesII)
     6671           
     6672    def test_read_mux_platform_problem3(self):
     6673       
     6674        # This is to test a situation where read_mux returned
     6675        # wrong values Win32
     6676
     6677       
     6678        from urs_ext import read_mux2
     6679       
     6680        from anuga.config import single_precision as epsilon       
     6681       
     6682        verbose = False
    66766683               
    6677                     # South is positive in MUX
     6684        tide = 1.5
     6685        time_step_count = 10
     6686        time_step = 0.02
     6687
     6688        '''
     6689        Win results
     6690        time_step = 0.2000001
     6691        This is OK       
     6692        '''
     6693       
     6694        '''
     6695        Win results
     6696        time_step = 0.20000001
     6697
     6698        ======================================================================
     6699ERROR: test_read_mux_platform_problem3 (__main__.Test_Data_Manager)
     6700----------------------------------------------------------------------
     6701Traceback (most recent call last):
     6702  File "test_data_manager.py", line 6718, in test_read_mux_platform_problem3
     6703    ha1[0]=num.sin(times_ref)
     6704ValueError: matrices are not aligned for copy
     6705
     6706        '''
     6707
     6708        '''
     6709        Win results
     6710        time_step = 0.200000001
     6711        FAIL
     6712         assert num.allclose(data[i][:parameters_index],
     6713         -va1[permutation[i], :])
     6714        '''
     6715        times_ref = num.arange(0, time_step_count*time_step, time_step)
     6716        #print "times_ref", times_ref
     6717       
     6718        lat_long_points = [(-21.5,114.5), (-21,114.5), (-21.5,115),
     6719                           (-21.,115.), (-22., 117.)]
     6720        stations = len(lat_long_points)
     6721       
     6722        # Create different timeseries starting and ending at different times
     6723        first_tstep=num.ones(stations, num.int)
     6724        first_tstep[0]+=2   # Point 0 starts at 2
     6725        first_tstep[1]+=4   # Point 1 starts at 4       
     6726        first_tstep[2]+=3   # Point 2 starts at 3
     6727       
     6728        last_tstep=(time_step_count)*num.ones(stations, num.int)
     6729        last_tstep[0]-=1    # Point 0 ends 1 step early
     6730        last_tstep[1]-=2    # Point 1 ends 2 steps early               
     6731        last_tstep[4]-=3    # Point 4 ends 3 steps early       
     6732       
     6733        # Create varying elevation data (positive values for seafloor)
     6734        gauge_depth=20*num.ones(stations, num.float)
     6735        for i in range(stations):
     6736            gauge_depth[i] += i**2
     6737           
     6738        # Create data to be written to second mux file       
     6739        ha1=num.ones((stations,time_step_count), num.float)
     6740        ha1[0]=num.sin(times_ref)
     6741        ha1[1]=2*num.sin(times_ref - 3)
     6742        ha1[2]=5*num.sin(4*times_ref)
     6743        ha1[3]=num.sin(times_ref)
     6744        ha1[4]=num.sin(2*times_ref-0.7)
     6745               
     6746        ua1=num.zeros((stations,time_step_count),num.float)
     6747        ua1[0]=3*num.cos(times_ref)       
     6748        ua1[1]=2*num.sin(times_ref-0.7)   
     6749        ua1[2]=num.arange(3*time_step_count,4*time_step_count)
     6750        ua1[4]=2*num.ones(time_step_count)
     6751       
     6752        va1=num.zeros((stations,time_step_count),num.float)
     6753        va1[0]=2*num.cos(times_ref-0.87)       
     6754        va1[1]=3*num.ones(time_step_count)
     6755        va1[3]=2*num.sin(times_ref-0.71)       
     6756        #print "va1[0]", va1[0]  # The 8th element is what will go bad.
     6757        # Ensure data used to write mux file to be zero when gauges are
     6758        # not recording
     6759        for i in range(stations):
     6760             # For each point
     6761             for j in range(0, first_tstep[i]-1) + range(last_tstep[i],
     6762                                                         time_step_count):
     6763                 # For timesteps before and after recording range
     6764                 ha1[i][j] = ua1[i][j] = va1[i][j] = 0.0
     6765
     6766
     6767        #print 'Second station to be written to MUX'
     6768        #print 'ha', ha1[0,:]
     6769        #print 'ua', ua1[0,:]
     6770        #print 'va', va1[0,:]
     6771       
     6772        # Write second mux file to be combined by urs2sts
     6773        base_nameII, filesII = self.write_mux2(lat_long_points,
     6774                                               time_step_count, time_step,
     6775                                               first_tstep, last_tstep,
     6776                                               depth=gauge_depth,
     6777                                               ha=ha1,
     6778                                               ua=ua1,
     6779                                               va=va1)
     6780        #print "filesII", filesII
     6781
     6782
     6783
     6784
     6785        # Read mux file back and verify it's correcness
     6786
     6787        ####################################################
     6788        # FIXME (Ole): This is where the test should
     6789        # verify that the MUX files are correct.
     6790
     6791        #JJ: It appears as though
     6792        #that certain quantities are not being stored with enough precision
     6793        #inn muxfile or more likely that they are being cast into a
     6794        #lower precision when read in using read_mux2 Time step and q_time
     6795        # are equal but only to approx 1e-7
     6796        ####################################################
     6797
     6798        #define information as it should be stored in mus2 files
     6799        points_num=len(lat_long_points)
     6800        depth=gauge_depth
     6801        ha=ha1
     6802        ua=ua1
     6803        va=va1
     6804       
     6805        quantities = ['HA','UA','VA']
     6806        mux_names = [WAVEHEIGHT_MUX2_LABEL,
     6807                     EAST_VELOCITY_MUX2_LABEL,
     6808                     NORTH_VELOCITY_MUX2_LABEL]
     6809        quantities_init = [[],[],[]]
     6810        latlondeps = []
     6811        #irrelevant header information
     6812        ig=ilon=ilat=0
     6813        mcolat=mcolon=centerlat=centerlon=offset=az=baz=id=0.0
     6814        # urs binary is latitude fastest
     6815        for i,point in enumerate(lat_long_points):
     6816            lat = point[0]
     6817            lon = point[1]
     6818            _ , e, n = redfearn(lat, lon)
     6819            if depth is None:
     6820                this_depth = n
     6821            else:
     6822                this_depth = depth[i]
     6823            latlondeps.append([lat, lon, this_depth])
     6824
     6825            if ha is None:
     6826                this_ha = e
     6827                quantities_init[0].append(num.ones(time_step_count,
     6828                                                   num.float)*this_ha) # HA
     6829            else:
     6830                quantities_init[0].append(ha[i])
     6831            if ua is None:
     6832                this_ua = n
     6833                quantities_init[1].append(num.ones(time_step_count,
     6834                                                   num.float)*this_ua) # UA
     6835            else:
     6836                quantities_init[1].append(ua[i])
     6837            if va is None:
     6838                this_va = e
     6839                quantities_init[2].append(num.ones(time_step_count,
     6840                                                   num.float)*this_va) #
     6841            else:
     6842                quantities_init[2].append(va[i])
     6843
     6844        for i, q in enumerate(quantities):
     6845            #print
     6846            #print i, q
     6847           
     6848            q_time = num.zeros((time_step_count, points_num), num.float)
     6849            quantities_init[i] = ensure_numeric(quantities_init[i])
     6850            for time in range(time_step_count):
     6851                #print i, q, time, quantities_init[i][:,time]
     6852                q_time[time,:] = quantities_init[i][:,time]
     6853                #print i, q, time, q_time[time, :]
     6854
     6855           
     6856            filename = base_nameII + mux_names[i]
     6857            f = open(filename, 'rb')
     6858            if self.verbose: print 'Reading' + filename
     6859            assert abs(points_num-unpack('i',f.read(4))[0])<epsilon
     6860            #write mux 2 header
     6861            for latlondep in latlondeps:
     6862                assert abs(latlondep[0]-unpack('f',f.read(4))[0])<epsilon
     6863                assert abs(latlondep[1]-unpack('f',f.read(4))[0])<epsilon
     6864                assert abs(mcolat-unpack('f',f.read(4))[0])<epsilon
     6865                assert abs(mcolon-unpack('f',f.read(4))[0])<epsilon
     6866                assert abs(ig-unpack('i',f.read(4))[0])<epsilon
     6867                assert abs(ilon-unpack('i',f.read(4))[0])<epsilon
     6868                assert abs(ilat-unpack('i',f.read(4))[0])<epsilon
     6869                assert abs(latlondep[2]-unpack('f',f.read(4))[0])<epsilon
     6870                assert abs(centerlat-unpack('f',f.read(4))[0])<epsilon
     6871                assert abs(centerlon-unpack('f',f.read(4))[0])<epsilon
     6872                assert abs(offset-unpack('f',f.read(4))[0])<epsilon
     6873                assert abs(az-unpack('f',f.read(4))[0])<epsilon
     6874                assert abs(baz-unpack('f',f.read(4))[0])<epsilon
     6875               
     6876                x = unpack('f', f.read(4))[0]
     6877                #print time_step
     6878                #print x
     6879                assert abs(time_step-x)<epsilon
     6880                assert abs(time_step_count-unpack('i',f.read(4))[0])<epsilon
     6881                for j in range(4): # identifier
     6882                    assert abs(id-unpack('i',f.read(4))[0])<epsilon
     6883
     6884            #first_tstep=1
     6885            #last_tstep=time_step_count
     6886            for i,latlondep in enumerate(latlondeps):
     6887                assert abs(first_tstep[i]-unpack('i',f.read(4))[0])<epsilon
     6888            for i,latlondep in enumerate(latlondeps):
     6889                assert abs(last_tstep[i]-unpack('i',f.read(4))[0])<epsilon
     6890
     6891            # Find when first station starts recording
     6892            min_tstep = min(first_tstep)
     6893            # Find when all stations have stopped recording
     6894            max_tstep = max(last_tstep)
     6895
     6896            #for time in  range(time_step_count):
     6897            for time in range(min_tstep-1,max_tstep):
     6898                assert abs(time*time_step-unpack('f',f.read(4))[0])<epsilon
     6899                for point_i in range(points_num):
     6900                    if time+1>=first_tstep[point_i] and time+1<=last_tstep[point_i]:
     6901                        x = unpack('f',f.read(4))[0]
     6902                        #print time, x, q_time[time, point_i]
     6903                        if q == 'VA': x = -x # South is positive in MUX
     6904                        #print q+" q_time[%d, %d] = %f" %(time, point_i,
     6905                                                      #q_time[time, point_i])
     6906                        assert abs(q_time[time, point_i]-x)<epsilon
     6907
     6908            f.close()
     6909                           
     6910        permutation = ensure_numeric([4,0,2])
     6911                   
     6912        # Create ordering file
     6913#         _, ordering_filename = tempfile.mkstemp('')
     6914#         order_fid = open(ordering_filename, 'w') 
     6915#         order_fid.write('index, longitude, latitude\n')
     6916#         for index in permutation:
     6917#             order_fid.write('%d, %f, %f\n' %(index,
     6918#                                              lat_long_points[index][1],
     6919#                                              lat_long_points[index][0]))
     6920#         order_fid.close()
     6921       
     6922        # -------------------------------------
     6923        # Now read files back and check values
     6924        weights = ensure_numeric([1.0])
     6925
     6926        # For each quantity read the associated list of source mux2 file with
     6927        # extention associated with that quantity
     6928        file_params=-1*num.ones(3,num.float) # [nsta,dt,nt]
     6929        OFFSET = 5
     6930
     6931        for j, file in enumerate(filesII):
     6932            # Read stage, u, v enumerated as j
     6933            #print 'Reading', j, file
     6934            #print "file", file
     6935            data = read_mux2(1, [file], weights, file_params,
     6936                             permutation, verbose)
     6937            #print str(j) + "data", data
     6938
     6939            #print 'Data received by Python'
     6940            #print data[1][8]
     6941            number_of_selected_stations = data.shape[0]
     6942            #print "number_of_selected_stations", number_of_selected_stations
     6943            #print "stations", stations
     6944
     6945            # Index where data ends and parameters begin
     6946            parameters_index = data.shape[1]-OFFSET         
     6947                 
     6948            for i in range(number_of_selected_stations):
     6949       
     6950                #print i, parameters_index
     6951                if j == 0:
     6952                    assert num.allclose(data[i][:parameters_index],
     6953                                        ha1[permutation[i], :])
     6954                   
     6955                if j == 1: assert num.allclose(data[i][:parameters_index], ua1[permutation[i], :])
     6956                if j == 2:
    66786957                    assert num.allclose(data[i][:parameters_index], -va1[permutation[i], :])
    6679                    
     6958       
     6959        self.delete_mux(filesII)           
    66806960       
    66816961    def test_urs2sts0(self):
     
    68197099
    68207100        base_name, files = self.write_mux2(lat_long_points,
    6821                                       time_step_count, time_step,
    6822                                       first_tstep, last_tstep,
    6823                                       depth=gauge_depth,
    6824                                       ha=ha,
    6825                                       ua=ua,
    6826                                       va=va)
     7101                                           time_step_count, time_step,
     7102                                           first_tstep, last_tstep,
     7103                                           depth=gauge_depth,
     7104                                           ha=ha,
     7105                                           ua=ua,
     7106                                           va=va)
    68277107
    68287108        urs2sts(base_name,
     
    68507130        y = points[:,1]
    68517131
    6852         #Check that all coordinate are correctly represented       
    6853         #Using the non standard projection (50)
     7132        # Check that all coordinate are correctly represented       
     7133        # Using the non standard projection (50)
    68547134        for i in range(4):
    68557135            zone, e, n = redfearn(lat_long_points[i][0],
     
    68587138            assert num.allclose([x[i],y[i]], [e,n])
    68597139            assert zone==-1
     7140       
     7141        self.delete_mux(files)
    68607142
    68617143           
     
    69157197        y = points[:,1]
    69167198
    6917         #Check that all coordinate are correctly represented       
    6918         #Using the non standard projection (50)
     7199        # Check that all coordinate are correctly represented       
     7200        # Using the non standard projection (50)
    69197201        for i in range(4):
    69207202            zone, e, n = redfearn(lat_long_points[i][0], lat_long_points[i][1],
     
    69227204            assert num.allclose([x[i],y[i]], [e,n])
    69237205            assert zone==geo_reference.zone
     7206       
     7207        self.delete_mux(files)
    69247208
    69257209           
     
    76527936            raise Exception, msg
    76537937
     7938       
     7939        self.delete_mux(filesI)
     7940        self.delete_mux(filesII)
    76547941
    76557942       
     
    80348321        os.remove(sts_file)
    80358322       
    8036        
    8037        
    80388323        #----------------------
    80398324        # Then read the mux files together and test
     
    81488433        os.remove(sts_file)
    81498434       
    8150 
    8151        
    81528435        #---------------
    81538436        # "Manually" add the timeseries up with weights and test
     
    81578440        stage_man = weights[0]*(stage0-tide) + weights[1]*(stage1-tide) + tide
    81588441        assert num.allclose(stage_man, stage)
    8159        
    8160        
    8161        
    8162        
    8163        
    81648442               
    81658443       
     
    82888566                       
    82898567        assert num.allclose(domain_fbound.quantities['xmomentum'].vertex_values,
    8290                             domain_drchlt.quantities['xmomentum'].vertex_values)                        
     8568                            domain_drchlt.quantities['xmomentum'].vertex_values)
    82918569                       
    82928570        assert num.allclose(domain_fbound.quantities['ymomentum'].vertex_values,
    8293                             domain_drchlt.quantities['ymomentum'].vertex_values)                                               
     8571                            domain_drchlt.quantities['ymomentum'].vertex_values)
    82948572       
    82958573       
    82968574        os.remove(sts_file+'.sts')
    82978575        os.remove(meshname)
    8298        
    8299        
    8300        
    83018576               
    83028577       
     
    84198694
    84208695
    8421 
    8422        
    8423        
    8424        
    8425            
    8426 
    8427        
    84288696        domain_drchlt = Domain(meshname)
    84298697        domain_drchlt.set_quantity('stage', tide)
     
    84478715                       
    84488716        assert num.allclose(domain_fbound.quantities['xmomentum'].vertex_values,
    8449                             domain_drchlt.quantities['xmomentum'].vertex_values)                        
     8717                            domain_drchlt.quantities['xmomentum'].vertex_values)
    84508718                       
    84518719        assert num.allclose(domain_fbound.quantities['ymomentum'].vertex_values,
    8452                             domain_drchlt.quantities['ymomentum'].vertex_values)                                               
    8453        
     8720                            domain_drchlt.quantities['ymomentum'].vertex_values)
    84548721       
    84558722        os.remove(sts_file+'.sts')
    84568723        os.remove(meshname)
    8457 
    8458 
    8459        
    84608724               
    84618725       
     
    1093411198        domain.set_quantity('xmomentum', uh)
    1093511199        domain.set_boundary( {'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})
     11200
    1093611201        for t in domain.evolve(yieldstep=1, finaltime = t_end):
    1093711202            pass
     
    1135711622
    1135811623           
     11624    def test_copy_code_files(self):
     11625        '''test that the copy_code_files() function is sane.'''
     11626
     11627        def create_file(f):
     11628            fd = open(f, 'w')
     11629            fd.write('%s\n' % f)
     11630            fd.close()
     11631
     11632        # create working directories and test files
     11633        work_dir = tempfile.mkdtemp()
     11634        dst_dir = tempfile.mkdtemp(dir=work_dir)
     11635        src_dir = tempfile.mkdtemp(dir=work_dir)
     11636
     11637        f1 = 'file1'       
     11638        filename1 = os.path.join(src_dir, f1)
     11639        create_file(filename1)
     11640        f2 = 'file2'       
     11641        filename2 = os.path.join(src_dir, f2)
     11642        create_file(filename2)
     11643        f3 = 'file3'       
     11644        filename3 = os.path.join(src_dir, f3)
     11645        create_file(filename3)
     11646        f4 = 'file4'       
     11647        filename4 = os.path.join(src_dir, f4)
     11648        create_file(filename4)
     11649        f5 = 'file5'       
     11650        filename5 = os.path.join(src_dir, f5)
     11651        create_file(filename5)
     11652
     11653        # exercise the copy function
     11654        copy_code_files(dst_dir, filename1)
     11655        copy_code_files(dst_dir, filename1, filename2)
     11656        copy_code_files(dst_dir, (filename4, filename5, filename3))
     11657
     11658        # test that files were actually copied
     11659        self.failUnless(access(os.path.join(dst_dir, f1), F_OK))
     11660        self.failUnless(access(os.path.join(dst_dir, f2), F_OK))
     11661        self.failUnless(access(os.path.join(dst_dir, f3), F_OK))
     11662        self.failUnless(access(os.path.join(dst_dir, f4), F_OK))
     11663        self.failUnless(access(os.path.join(dst_dir, f5), F_OK))
     11664
     11665        # clean up
     11666        shutil.rmtree(work_dir)
    1135911667           
    1136011668#-------------------------------------------------------------
     
    1136211670if __name__ == "__main__":
    1136311671    suite = unittest.makeSuite(Test_Data_Manager,'test')
    11364 
    11365     # FIXME (Ole): This is the test that fails under Windows
    11366     #suite = unittest.makeSuite(Test_Data_Manager,'test_read_mux_platform_problem2')
    11367     #suite = unittest.makeSuite(Test_Data_Manager,'test_file_boundary_stsIV')
    1136811672   
    1136911673    if len(sys.argv) > 1 and sys.argv[1][0].upper() == 'V':
  • branches/numpy/anuga/shallow_water/test_shallow_water_domain.py

    r6689 r6902  
    2727def linear_function(point):
    2828    point = num.array(point)
    29     return point[:,0] + point[:,1]
    30 
     29    return point[:,0]+point[:,1]
    3130
    3231class Weir:
    3332    """Set a bathymetry for weir with a hole and a downstream gutter
    34 
    3533    x,y are assumed to be in the unit square
    3634    """
     
    4543        z = num.zeros(N, num.float)
    4644        for i in range(N):
    47             z[i] = -x[i] / 2    # General slope
    48 
    49             # Flattish bit to the left
     45            z[i] = -x[i]/2  #General slope
     46
     47            #Flattish bit to the left
    5048            if x[i] < 0.3:
    5149                z[i] = -x[i]/10
    5250
    53             # Weir
     51            #Weir
    5452            if x[i] >= 0.3 and x[i] < 0.4:
    5553                z[i] = -x[i]+0.9
    5654
    57             # Dip
     55            #Dip
    5856            x0 = 0.6
    5957            depth = -1.0
     
    6260                if x[i] > x0 and x[i] < 0.9:
    6361                    z[i] = depth
    64                 # RHS plateaux
     62                #RHS plateaux
    6563                if x[i] >= 0.9:
    6664                    z[i] = plateaux
    6765            elif y[i] >= 0.7 and y[i] < 1.5:
    68                 # Restrict and deepen
     66                #Restrict and deepen
    6967                if x[i] >= x0 and x[i] < 0.8:
    7068                    z[i] = depth - (y[i]/3 - 0.3)
    7169                elif x[i] >= 0.8:
    72                     # RHS plateaux
     70                    #RHS plateaux
    7371                    z[i] = plateaux
    7472            elif y[i] >= 1.5:
    7573                if x[i] >= x0 and x[i] < 0.8 + (y[i]-1.5)/1.2:
    76                     # Widen up and stay at constant depth
    77                     z[i] = depth - 1.5/5
    78                 elif x[i] >= 0.8 + (y[i] - 1.5)/1.2:
    79                     # RHS plateaux
     74                    #Widen up and stay at constant depth
     75                    z[i] = depth-1.5/5
     76                elif x[i] >= 0.8 + (y[i]-1.5)/1.2:
     77                    #RHS plateaux
    8078                    z[i] = plateaux
    8179
    82             # Hole in weir (slightly higher than inflow condition)
     80            #Hole in weir (slightly higher than inflow condition)
    8381            if x[i] >= 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4:
    84                 z[i] = -x[i] + self.inflow_stage + 0.02
    85 
    86             # Channel behind weir
     82                z[i] = -x[i]+self.inflow_stage + 0.02
     83
     84            #Channel behind weir
    8785            x0 = 0.5
    8886            if x[i] >= 0.4 and x[i] < x0 and y[i] > 0.2 and y[i] < 0.4:
    89                 z[i] = -x[i] + self.inflow_stage + 0.02
     87                z[i] = -x[i]+self.inflow_stage + 0.02
    9088
    9189            if x[i] >= x0 and x[i] < 0.6 and y[i] > 0.2 and y[i] < 0.4:
    92                 # Flatten it out towards the end
    93                 z[i] = -x0 + self.inflow_stage + 0.02 + (x0 - x[i])/5
     90                #Flatten it out towards the end
     91                z[i] = -x0+self.inflow_stage + 0.02 + (x0-x[i])/5
    9492
    9593            # Hole to the east
     
    9795            y0 = 0.35
    9896            if num.sqrt((2*(x[i]-x0))**2 + (2*(y[i]-y0))**2) < 0.2:
    99                 z[i] = num.sqrt((x[i]-x0)**2 + (y[i]-y0)**2) - 1.0
    100 
    101             # Tiny channel draining hole
     97                z[i] = num.sqrt(((x[i]-x0))**2 + ((y[i]-y0))**2)-1.0
     98
     99            #Tiny channel draining hole
    102100            if x[i] >= 1.14 and x[i] < 1.2 and y[i] >= 0.4 and y[i] < 0.6:
    103                 z[i] = -0.9    # North south
     101                z[i] = -0.9 #North south
    104102
    105103            if x[i] >= 0.9 and x[i] < 1.18 and y[i] >= 0.58 and y[i] < 0.65:
    106                 z[i] = -1.0 + (x[i]-0.9)/3    # East west
     104                z[i] = -1.0 + (x[i]-0.9)/3 #East west
    107105
    108106            # Stuff not in use
     
    136134        z = num.zeros(N, num.float)
    137135        for i in range(N):
    138             z[i] = -x[i]    # General slope
    139 
    140             # Flat bit to the left
     136            z[i] = -x[i]  #General slope
     137
     138            #Flat bit to the left
    141139            if x[i] < 0.3:
    142                 z[i] = -x[i]/10    # General slope
    143 
    144             # Weir
     140                z[i] = -x[i]/10  #General slope
     141
     142            #Weir
    145143            if x[i] > 0.3 and x[i] < 0.4:
    146                 z[i] = -x[i] + 0.9
    147 
    148             # Dip
     144                z[i] = -x[i]+0.9
     145
     146            #Dip
    149147            if x[i] > 0.6 and x[i] < 0.9:
    150                 z[i] = -x[i] - 0.5    # -y[i]/5
    151 
    152             # Hole in weir (slightly higher than inflow condition)
     148                z[i] = -x[i]-0.5  #-y[i]/5
     149
     150            #Hole in weir (slightly higher than inflow condition)
    153151            if x[i] > 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4:
    154                 z[i] = -x[i] + self.inflow_stage + 0.05
     152                z[i] = -x[i]+self.inflow_stage + 0.05
     153
    155154
    156155        return z/2
     
    170169
    171170    N = len(x)
    172     s = 0 * x    # New array
     171    s = 0*x  #New array
    173172
    174173    for k in range(N):
     
    263262        H0 = 0.0
    264263
    265         max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
    266                                   epsilon, g, H0)
    267 
    268         assert num.allclose(edgeflux, [0, 0, 0])
     264       
     265        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0)
     266
     267        assert num.allclose(edgeflux, [0,0,0])
    269268        assert max_speed == 0.
    270269
     
    272271        w = 2.0
    273272
    274         normal = num.array([1., 0])
     273        normal = num.array([1.,0])
    275274        ql = num.array([w, 0, 0])
    276275        qr = num.array([w, 0, 0])
     
    280279        H0 = 0.0
    281280
    282         max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
    283                                   epsilon, g, H0)
    284 
     281        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0)       
    285282        assert num.allclose(edgeflux, [0., 0.5*g*h**2, 0.])
    286283        assert max_speed == num.sqrt(g*h)
     
    290287    #    w = 2.0
    291288    #
    292     #    normal = array([1., 0])
     289    #    normal = array([1.,0])
    293290    #    ql = array([w, 0, 0])
    294291    #    qr = array([w, 0, 0])
     
    15171514        uh = u*h
    15181515
    1519         Br = Reflective_boundary(domain)       # Side walls
    1520         Bd = Dirichlet_boundary([w, uh, 0])    # 2 m/s across the 3 m inlet:
     1516        Br = Reflective_boundary(domain)     # Side walls
     1517        Bd = Dirichlet_boundary([w, uh, 0])  # 2 m/s across the 3 m inlet:
     1518
    15211519
    15221520        # Initial conditions
  • branches/numpy/anuga/shallow_water/urs_ext.c

    r6780 r6902  
    104104                memcpy(data + it, muxData + offset, sizeof(float));
    105105               
    106                 //printf("%d: muxdata=%f\n", it, muxData[offset]);                     
    107                 //printf("data[%d]=%f, offset=%d\n", it, data[it], offset);       
     106                //printf("%d: muxdata=%f\n", it, muxData[offset]); 
     107                //printf("data[%d]=%f, offset=%d\n", it, data[it], offset);
    108108                offset++;
    109109            }
     
    221221
    222222        // Open the mux file
    223         if((fp = fopen(muxFileName, "r")) == NULL)
     223        if((fp = fopen(muxFileName, "rb")) == NULL)
    224224        {
    225225            char *err_msg = strerror(errno);
     
    237237            }
    238238       
    239             fros = (int*) malloc(*total_number_of_stations*numSrc*sizeof(int)); 
     239            fros = (int*) malloc(*total_number_of_stations*numSrc*sizeof(int));
    240240            lros = (int*) malloc(*total_number_of_stations*numSrc*sizeof(int));
    241241     
     
    447447    temp_sts_data = (float*) calloc(len_sts_data, sizeof(float));
    448448
    449     muxData = (float*) malloc(numDataMax*sizeof(float));
     449    muxData = (float*) calloc(numDataMax, sizeof(float));
    450450   
    451451    // Loop over all sources
     
    460460        // Read in data block from mux2 file
    461461        muxFileName = muxFileNameArray[isrc];
    462         if((fp = fopen(muxFileName, "r")) == NULL)
     462        if((fp = fopen(muxFileName, "rb")) == NULL)
    463463        {
    464464            fprintf(stderr, "cannot open file %s\n", muxFileName);
     
    470470        }
    471471
    472         offset = sizeof(int) + total_number_of_stations*(sizeof(struct tgsrwg) + 2*sizeof(int));
    473         fseek(fp, offset, 0);
     472        offset = (long int)sizeof(int) + total_number_of_stations*(sizeof(struct tgsrwg) + 2*sizeof(int));
     473        //printf("\n offset %i ", (long int)offset);
     474        fseek(fp, offset, 0);
    474475
    475476        numData = getNumData(fros_per_source,
    476477                             lros_per_source,
    477478                             total_number_of_stations);
    478                              
    479                              
    480                                      
    481         elements_read = fread(muxData, ((int) numData)*sizeof(float), 1, fp);
     479        // Note numData is larger than what it has to be.                   
     480        //elements_read = fread(muxData, ((int) numData)*sizeof(float), 1, fp);
     481        elements_read = fread(muxData, (size_t) sizeof(float), (size_t) numData, fp);
     482        //printf("\n elements_read  %d, ", (int)elements_read);
     483        //printf("\n ferror(fp)  %d, ", (int)ferror(fp));
    482484        if ((int) elements_read == 0 && ferror(fp)) {
    483485          fprintf(stderr, "Error reading mux data\n");
    484486          return NULL;
    485         }
    486                
    487         fclose(fp);
     487        }       
    488488       
    489 
    490         // FIXME (Ole): This is where Nariman and Ole traced the platform dependent
    491         // difference on 11 November 2008. We don't think the problem lies in the
    492         // C code. Maybe it is a problem with the MUX files written by the unit test
    493         // that fails on Windows but works OK on Linux. JJ's test on 17th November shows
    494         // that as far as Python is concerned this file should be OK on both platforms.
    495        
    496        
    497        
    498         // These printouts are enough to show the problem when compared
    499         // on the two platforms
    500         //printf("\nRead %d elements, ", (int) numData);
    501         //printf("muxdata[%d]=%f\n", 39, muxData[39]);         
    502        
    503        
    504         /*
    505         In Windows we get
    506        
    507         Read 85 elements, muxdata[39]=0.999574
    508         Read 85 elements, muxdata[39]=-0.087599
    509         Read 85 elements, muxdata[39]=-0.087599
    510        
    511        
    512         In Linux we get (the correct)
    513        
    514         Read 85 elements, muxdata[39]=0.999574
    515         Read 85 elements, muxdata[39]=-0.087599
    516         Read 85 elements, muxdata[39]=1.490349
    517         */
    518        
    519        
    520        
    521        
     489        fclose(fp); 
    522490       
    523491        // loop over stations present in the permutation array
Note: See TracChangeset for help on using the changeset viewer.