Changeset 6737


Ignore:
Timestamp:
Apr 7, 2009, 9:38:40 AM (16 years ago)
Author:
nariman
Message:

Refactoring of code in _flux_function_central for optimisation.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r6703 r6737  
    5252
    5353int find_qmin_and_qmax(double dq0, double dq1, double dq2,
    54                        double *qmin, double *qmax){
     54               double *qmin, double *qmax){
    5555  // Considering the centroid of an FV triangle and the vertices of its
    5656  // auxiliary triangle, find
    … …  
    6767    if (dq1>=dq2){
    6868      if (dq1>=0.0)
    69         *qmax=dq0+dq1;
     69    *qmax=dq0+dq1;
    7070      else
    71         *qmax=dq0;
    72        
     71    *qmax=dq0;
     72   
    7373      *qmin=dq0+dq2;
    7474      if (*qmin>=0.0) *qmin = 0.0;
    … …  
    7676    else{// dq1<dq2
    7777      if (dq2>0)
    78         *qmax=dq0+dq2;
     78    *qmax=dq0+dq2;
    7979      else
    80         *qmax=dq0;
    81        
    82       *qmin=dq0+dq1;   
     80    *qmax=dq0;
     81   
     82      *qmin=dq0+dq1;   
    8383      if (*qmin>=0.0) *qmin=0.0;
    8484    }
    … …  
    8787    if (dq1<=dq2){
    8888      if (dq1<0.0)
    89         *qmin=dq0+dq1;
     89    *qmin=dq0+dq1;
    9090      else
    91         *qmin=dq0;
    92        
    93       *qmax=dq0+dq2;   
     91    *qmin=dq0;
     92   
     93      *qmax=dq0+dq2;   
    9494      if (*qmax<=0.0) *qmax=0.0;
    9595    }
    9696    else{// dq1>dq2
    9797      if (dq2<0.0)
    98         *qmin=dq0+dq2;
     98    *qmin=dq0+dq2;
    9999      else
    100         *qmin=dq0;
    101        
     100    *qmin=dq0;
     101   
    102102      *qmax=dq0+dq1;
    103103      if (*qmax<=0.0) *qmax=0.0;
    … …  
    141141
    142142double compute_froude_number(double uh,
    143                              double h,
    144                              double g,
    145                              double epsilon) {
    146                          
     143                 double h,
     144                 double g,
     145                 double epsilon) {
     146             
    147147  // Compute Froude number; v/sqrt(gh)
    148148  // FIXME (Ole): Not currently in use
    … …  
    169169// FIXME: Perhaps inline this and profile
    170170double _compute_speed(double *uh,
    171                       double *h,
    172                       double epsilon,
    173                       double h0) {
     171              double *h,
     172              double epsilon,
     173              double h0) {
    174174 
    175175  double u;
    … …  
    192192// Innermost flux function (using stage w=z+h)
    193193int _flux_function_central(double *q_left, double *q_right,
    194                            double z_left, double z_right,
    195                            double n1, double n2,
    196                            double epsilon, double H0, double g,
    197                            double *edgeflux, double *max_speed) {
     194                           double z_left, double z_right,
     195                           double n1, double n2,
     196                           double epsilon, double H0, double g,
     197                           double *edgeflux, double *max_speed)
     198{
    198199
    199200  /*Compute fluxes between volumes for the shallow water wave equation
    … …  
    221222  double h0 = H0*H0; // This ensures a good balance when h approaches H0.
    222223                     // But evidence suggests that h0 can be as little as
    223                      // epsilon!
     224             // epsilon!
    224225 
    225226  // Copy conserved quantities to protect from modification
    226   for (i=0; i<3; i++) {
    227     q_left_rotated[i] = q_left[i];
    228     q_right_rotated[i] = q_right[i];
    229   }
     227  q_left_rotated[0] = q_left[0];
     228  q_right_rotated[0] = q_right[0];
     229  q_left_rotated[1] = q_left[1];
     230  q_right_rotated[1] = q_right[1];
     231  q_left_rotated[2] = q_left[2];
     232  q_right_rotated[2] = q_right[2];
    230233
    231234  // Align x- and y-momentum with x-axis
    … …  
    233236  _rotate(q_right_rotated, n1, n2);
    234237
    235   z = 0.5*(z_left+z_right); // Average elevation values.
     238  z = 0.5*(z_left + z_right); // Average elevation values.
    236239                            // Even though this will nominally allow for discontinuities
    237                             // in the elevation data, there is currently no numerical
    238                             // support for this so results may be strange near jumps in the bed.
     240                            // in the elevation data, there is currently no numerical
     241                            // support for this so results may be strange near jumps in the bed.
    239242
    240243  // Compute speeds in x-direction
    241244  w_left = q_left_rotated[0];         
    242   h_left = w_left-z;
     245  h_left = w_left - z;
    243246  uh_left = q_left_rotated[1];
    244247  u_left = _compute_speed(&uh_left, &h_left, epsilon, h0);
    245248
    246249  w_right = q_right_rotated[0];
    247   h_right = w_right-z;
     250  h_right = w_right - z;
    248251  uh_right = q_right_rotated[1];
    249252  u_right = _compute_speed(&uh_right, &h_right, epsilon, h0); 
    … …  
    261264  soundspeed_right = sqrt(g*h_right);
    262265
    263   s_max = max(u_left+soundspeed_left, u_right+soundspeed_right);
    264   if (s_max < 0.0) s_max = 0.0;
    265 
    266   s_min = min(u_left-soundspeed_left, u_right-soundspeed_right);
    267   if (s_min > 0.0) s_min = 0.0;
    268 
     266  s_max = max(u_left + soundspeed_left, u_right + soundspeed_right);
     267  if (s_max < 0.0)
     268  {
     269    s_max = 0.0;
     270  }
     271
     272  s_min = min(u_left - soundspeed_left, u_right - soundspeed_right);
     273  if (s_min > 0.0)
     274  {
     275    s_min = 0.0;
     276  }
    269277 
    270278  // Flux formulas
    … …  
    277285  flux_right[2] = u_right*vh_right;
    278286
    279  
    280287  // Flux computation
    281   denom = s_max-s_min;
    282   if (denom < epsilon) { // FIXME (Ole): Try using H0 here
    283     for (i=0; i<3; i++) edgeflux[i] = 0.0;
     288  denom = s_max - s_min;
     289  if (denom < epsilon)
     290  { // FIXME (Ole): Try using H0 here
     291    memset(edgeflux, 0, 3*sizeof(double));
    284292    *max_speed = 0.0;
    285   } else {
     293  }
     294  else
     295  {
    286296    inverse_denominator = 1.0/denom;
    287     for (i=0; i<3; i++) {
     297    for (i = 0; i < 3; i++)
     298    {
    288299      edgeflux[i] = s_max*flux_left[i] - s_min*flux_right[i];
    289       edgeflux[i] += s_max*s_min*(q_right_rotated[i]-q_left_rotated[i]);
     300      edgeflux[i] += s_max*s_min*(q_right_rotated[i] - q_left_rotated[i]);
    290301      edgeflux[i] *= inverse_denominator;
    291302    }
    … …  
    318329// Computational function for flux computation (using stage w=z+h)
    319330int flux_function_kinetic(double *q_left, double *q_right,
    320                   double z_left, double z_right,
    321                   double n1, double n2,
    322                   double epsilon, double H0, double g,
    323                   double *edgeflux, double *max_speed) {
     331          double z_left, double z_right,
     332          double n1, double n2,
     333          double epsilon, double H0, double g,
     334          double *edgeflux, double *max_speed) {
    324335
    325336  /*Compute fluxes between volumes for the shallow water wave equation
    … …  
    416427
    417428void _manning_friction(double g, double eps, int N,
    418                        double* w, double* z,
    419                        double* uh, double* vh,
    420                        double* eta, double* xmom, double* ymom) {
     429               double* w, double* z,
     430               double* uh, double* vh,
     431               double* eta, double* xmom, double* ymom) {
    421432
    422433  int k;
    … …  
    444455/*
    445456void _manning_friction_explicit(double g, double eps, int N,
    446                        double* w, double* z,
    447                        double* uh, double* vh,
    448                        double* eta, double* xmom, double* ymom) {
     457               double* w, double* z,
     458               double* uh, double* vh,
     459               double* eta, double* xmom, double* ymom) {
    449460
    450461  int k;
    … …  
    455466      h = w[k]-z[k];
    456467      if (h >= eps) {
    457         S = -g * eta[k]*eta[k] * sqrt((uh[k]*uh[k] + vh[k]*vh[k]));
    458         S /= pow(h, 7.0/3);      //Expensive (on Ole's home computer)
    459         //S /= exp(7.0/3.0*log(h));      //seems to save about 15% over manning_friction
    460         //S /= h*h*(1 + h/3.0 - h*h/9.0); //FIXME: Could use a Taylor expansion
    461 
    462 
    463         //Update momentum
    464         xmom[k] += S*uh[k];
    465         ymom[k] += S*vh[k];
     468    S = -g * eta[k]*eta[k] * sqrt((uh[k]*uh[k] + vh[k]*vh[k]));
     469    S /= pow(h, 7.0/3);      //Expensive (on Ole's home computer)
     470    //S /= exp(7.0/3.0*log(h));      //seems to save about 15% over manning_friction
     471    //S /= h*h*(1 + h/3.0 - h*h/9.0); //FIXME: Could use a Taylor expansion
     472
     473
     474    //Update momentum
     475    xmom[k] += S*uh[k];
     476    ymom[k] += S*vh[k];
    466477      }
    467478    }
    … …  
    474485
    475486double velocity_balance(double uh_i,
    476                         double uh,
    477                         double h_i,
    478                         double h,
    479                         double alpha,
    480                         double epsilon) {
     487            double uh,
     488            double h_i,
     489            double h,
     490            double alpha,
     491            double epsilon) {
    481492  // Find alpha such that speed at the vertex is within one
    482   // order of magnitude of the centroid speed   
     493  // order of magnitude of the centroid speed   
    483494
    484495  // FIXME(Ole): Work in progress
    … …  
    489500 
    490501  printf("alpha = %f, uh_i=%f, uh=%f, h_i=%f, h=%f\n",
    491         alpha, uh_i, uh, h_i, h);
     502    alpha, uh_i, uh, h_i, h);
    492503     
    493504 
    … …  
    503514  if (h < epsilon) {
    504515    b = 1.0e10; // Limit
    505   } else {     
     516  } else { 
    506517    b = m*fabs(h_i - h)/h;
    507518  }
    … …  
    510521
    511522  if (a > b) {
    512     estimate = (m-1)/(a-b);             
     523    estimate = (m-1)/(a-b);     
    513524   
    514525    printf("Alpha %f, estimate %f\n",
    515            alpha, estimate);   
    516            
     526       alpha, estimate);   
     527       
    517528    if (alpha < estimate) {
    518529      printf("Adjusting alpha from %f to %f\n",
    519              alpha, estimate);
     530         alpha, estimate);
    520531      alpha = estimate;
    521532    }
    … …  
    523534 
    524535    if (h < h_i) {
    525       estimate = (m-1)/(a-b);                
     536      estimate = (m-1)/(a-b);            
    526537   
    527538      printf("Alpha %f, estimate %f\n",
    528              alpha, estimate);   
    529            
     539         alpha, estimate);   
     540       
    530541      if (alpha < estimate) {
    531         printf("Adjusting alpha from %f to %f\n",
    532                alpha, estimate);
    533         alpha = estimate;
     542    printf("Adjusting alpha from %f to %f\n",
     543           alpha, estimate);
     544    alpha = estimate;
    534545      }   
    535546    }
    … …  
    543554
    544555int _balance_deep_and_shallow(int N,
    545                               double* wc,
    546                               double* zc,
    547                               double* wv,
    548                               double* zv,
    549                               double* hvbar, // Retire this
    550                               double* xmomc,
    551                               double* ymomc,
    552                               double* xmomv,
    553                               double* ymomv,
    554                               double H0,
    555                               int tight_slope_limiters,
    556                               int use_centroid_velocities,
    557                               double alpha_balance) {
     556                  double* wc,
     557                  double* zc,
     558                  double* wv,
     559                  double* zv,
     560                  double* hvbar, // Retire this
     561                  double* xmomc,
     562                  double* ymomc,
     563                  double* xmomv,
     564                  double* ymomv,
     565                  double H0,
     566                  int tight_slope_limiters,
     567                  int use_centroid_velocities,
     568                  double alpha_balance) {
    558569
    559570  int k, k3, i;
    … …  
    582593      // FIXME: Try with this one precomputed
    583594      for (i=0; i<3; i++) {
    584         dz = max(dz, fabs(zv[k3+i]-zc[k]));
     595    dz = max(dz, fabs(zv[k3+i]-zc[k]));
    585596      }
    586597    }
    … …  
    595606
    596607    //if (hmin < 0.0 ) {
    597     //  printf("hmin = %f\n",hmin);
     608    //  printf("hmin = %f\n",hmin);
    598609    //}
    599610   
    … …  
    610621     
    611622      if (dz > 0.0) {
    612         alpha = max( min( alpha_balance*hmin/dz, 1.0), 0.0 );     
     623    alpha = max( min( alpha_balance*hmin/dz, 1.0), 0.0 );     
    613624      } else {
    614         alpha = 1.0;  // Flat bed
     625    alpha = 1.0;  // Flat bed
    615626      }
    616627      //printf("Using old style limiter\n");
    … …  
    625636   
    626637      if (hmin < H0) {
    627         alpha = 1.0;
    628         for (i=0; i<3; i++) {
    629          
    630           h_diff = hc_k - hv[i];         
    631           if (h_diff <= 0) {
    632             // Deep water triangle is further away from bed than
    633             // shallow water (hbar < h). Any alpha will do
    634          
    635           } else { 
    636             // Denominator is positive which means that we need some of the
    637             // h-limited stage.
    638            
    639             alpha = min(alpha, (hc_k - H0)/h_diff);        
    640           }
    641         }
    642 
    643         // Ensure alpha in [0,1]
    644         if (alpha>1.0) alpha=1.0;
    645         if (alpha<0.0) alpha=0.0;
    646        
     638    alpha = 1.0;
     639    for (i=0; i<3; i++) {
     640     
     641      h_diff = hc_k - hv[i];     
     642      if (h_diff <= 0) {
     643        // Deep water triangle is further away from bed than
     644        // shallow water (hbar < h). Any alpha will do
     645     
     646      } else { 
     647        // Denominator is positive which means that we need some of the
     648        // h-limited stage.
     649       
     650        alpha = min(alpha, (hc_k - H0)/h_diff);    
     651      }
     652    }
     653
     654    // Ensure alpha in [0,1]
     655    if (alpha>1.0) alpha=1.0;
     656    if (alpha<0.0) alpha=0.0;
     657   
    647658      } else {
    648         // Use w-limited stage exclusively in deeper water.
    649         alpha = 1.0;       
     659    // Use w-limited stage exclusively in deeper water.
     660    alpha = 1.0;       
    650661      }
    651662    }
    652663
    653664
    654     //  Let
     665    //  Let
    655666    //
    656     //    wvi be the w-limited stage (wvi = zvi + hvi)
    657     //    wvi- be the h-limited state (wvi- = zvi + hvi-)
     667    //    wvi be the w-limited stage (wvi = zvi + hvi)
     668    //    wvi- be the h-limited state (wvi- = zvi + hvi-)
    658669    //
    659670    //
    660     //  where i=0,1,2 denotes the vertex ids
     671    //  where i=0,1,2 denotes the vertex ids
    661672    //
    662673    //  Weighted balance between w-limited and h-limited stage is
    663674    //
    664     //    wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi)
     675    //    wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi)
    665676    //
    666677    //  It follows that the updated wvi is
    667678    //    wvi := zvi + (1-alpha)*hvi- + alpha*hvi
    668679    //
    669     //  Momentum is balanced between constant and limited
     680    //  Momentum is balanced between constant and limited
    670681
    671682
    … …  
    673684      for (i=0; i<3; i++) {
    674685     
    675         wv[k3+i] = zv[k3+i] + (1-alpha)*hc_k + alpha*hv[i];     
    676 
    677         // Update momentum at vertices
    678         if (use_centroid_velocities == 1) {
    679           // This is a simple, efficient and robust option
    680           // It uses first order approximation of velocities, but retains
    681           // the order used by stage.
    682        
    683           // Speeds at centroids
    684           if (hc_k > epsilon) {
    685             uc = xmomc[k]/hc_k;
    686             vc = ymomc[k]/hc_k;
    687           } else {
    688             uc = 0.0;
    689             vc = 0.0;
    690           }
    691          
    692           // Vertex momenta guaranteed to be consistent with depth guaranteeing
    693           // controlled speed
    694           hv[i] = wv[k3+i] - zv[k3+i]; // Recompute (balanced) vertex depth
    695           xmomv[k3+i] = uc*hv[i];
    696           ymomv[k3+i] = vc*hv[i];
    697          
    698         } else {
    699           // Update momentum as a linear combination of
    700           // xmomc and ymomc (shallow) and momentum
    701           // from extrapolator xmomv and ymomv (deep).
    702           // This assumes that values from xmomv and ymomv have
    703           // been established e.g. by the gradient limiter.
    704 
    705           // FIXME (Ole): I think this should be used with vertex momenta
    706           // computed above using centroid_velocities instead of xmomc
    707           // and ymomc as they'll be more representative first order
    708           // values.
    709          
    710           xmomv[k3+i] = (1-alpha)*xmomc[k] + alpha*xmomv[k3+i];
    711           ymomv[k3+i] = (1-alpha)*ymomc[k] + alpha*ymomv[k3+i];
    712        
    713         }
     686    wv[k3+i] = zv[k3+i] + (1-alpha)*hc_k + alpha*hv[i];
     687
     688    // Update momentum at vertices
     689    if (use_centroid_velocities == 1) {
     690      // This is a simple, efficient and robust option
     691      // It uses first order approximation of velocities, but retains
     692      // the order used by stage.
     693   
     694      // Speeds at centroids
     695      if (hc_k > epsilon) {
     696        uc = xmomc[k]/hc_k;
     697        vc = ymomc[k]/hc_k;
     698      } else {
     699        uc = 0.0;
     700        vc = 0.0;
     701      }
     702     
     703      // Vertex momenta guaranteed to be consistent with depth guaranteeing
     704      // controlled speed
     705      hv[i] = wv[k3+i] - zv[k3+i]; // Recompute (balanced) vertex depth
     706      xmomv[k3+i] = uc*hv[i];
     707      ymomv[k3+i] = vc*hv[i];
     708     
     709    } else {
     710      // Update momentum as a linear combination of
     711      // xmomc and ymomc (shallow) and momentum
     712      // from extrapolator xmomv and ymomv (deep).
     713      // This assumes that values from xmomv and ymomv have
     714      // been established e.g. by the gradient limiter.
     715
     716      // FIXME (Ole): I think this should be used with vertex momenta
     717      // computed above using centroid_velocities instead of xmomc
     718      // and ymomc as they'll be more representative first order
     719      // values.
     720     
     721      xmomv[k3+i] = (1-alpha)*xmomc[k] + alpha*xmomv[k3+i];
     722      ymomv[k3+i] = (1-alpha)*ymomc[k] + alpha*ymomv[k3+i];
     723   
     724    }
    714725      }
    715726    }
    … …  
    722733
    723734int _protect(int N,
    724              double minimum_allowed_height,
    725              double maximum_allowed_speed,
    726              double epsilon,
    727              double* wc,
    728              double* zc,
    729              double* xmomc,
    730              double* ymomc) {
     735         double minimum_allowed_height,
     736         double maximum_allowed_speed,
     737         double epsilon,
     738         double* wc,
     739         double* zc,
     740         double* xmomc,
     741         double* ymomc) {
    731742
    732743  int k;
    … …  
    741752
    742753      if (hc < minimum_allowed_height) {
    743        
    744         // Set momentum to zero and ensure h is non negative
    745         xmomc[k] = 0.0;
    746         ymomc[k] = 0.0;
    747         if (hc <= 0.0) wc[k] = zc[k];
     754       
     755    // Set momentum to zero and ensure h is non negative
     756    xmomc[k] = 0.0;
     757    ymomc[k] = 0.0;
     758    if (hc <= 0.0) wc[k] = zc[k];
    748759      }
    749760    }
    … …  
    760771             
    761772        if (hc <= 0.0) {
    762                 wc[k] = zc[k];
    763         xmomc[k] = 0.0;
    764         ymomc[k] = 0.0;
     773            wc[k] = zc[k];
     774        xmomc[k] = 0.0;
     775        ymomc[k] = 0.0;
    765776        } else {
    766777          //Reduce excessive speeds derived from division by small hc
    767         //FIXME (Ole): This may be unnecessary with new slope limiters
    768         //in effect.
     778        //FIXME (Ole): This may be unnecessary with new slope limiters
     779        //in effect.
    769780         
    770781          u = xmomc[k]/hc;
    771           if (fabs(u) > maximum_allowed_speed) {
    772             reduced_speed = maximum_allowed_speed * u/fabs(u);
    773             //printf("Speed (u) has been reduced from %.3f to %.3f\n",
    774             //  u, reduced_speed);
    775             xmomc[k] = reduced_speed * hc;
    776           }
     782      if (fabs(u) > maximum_allowed_speed) {
     783        reduced_speed = maximum_allowed_speed * u/fabs(u);
     784        //printf("Speed (u) has been reduced from %.3f to %.3f\n",
     785        //  u, reduced_speed);
     786        xmomc[k] = reduced_speed * hc;
     787      }
    777788
    778789          v = ymomc[k]/hc;
    779           if (fabs(v) > maximum_allowed_speed) {
    780             reduced_speed = maximum_allowed_speed * v/fabs(v);
    781             //printf("Speed (v) has been reduced from %.3f to %.3f\n",
    782             //  v, reduced_speed);
    783             ymomc[k] = reduced_speed * hc;
    784           }
     790      if (fabs(v) > maximum_allowed_speed) {
     791        reduced_speed = maximum_allowed_speed * v/fabs(v);
     792        //printf("Speed (v) has been reduced from %.3f to %.3f\n",
     793        //  v, reduced_speed);
     794        ymomc[k] = reduced_speed * hc;
     795      }
    785796        }
    786797      }
    … …  
    793804
    794805int _assign_wind_field_values(int N,
    795                               double* xmom_update,
    796                               double* ymom_update,
    797                               double* s_vec,
    798                               double* phi_vec,
    799                               double cw) {
     806                  double* xmom_update,
     807                  double* ymom_update,
     808                  double* s_vec,
     809                  double* phi_vec,
     810                  double cw) {
    800811
    801812  // Assign windfield values to momentum updates
    … …  
    842853
    843854  if (!PyArg_ParseTuple(args, "OOOddOddd",
    844                         &normal, &ql, &qr, &zl, &zr, &edgeflux,
    845                         &epsilon, &g, &H0)) {
     855            &normal, &ql, &qr, &zl, &zr, &edgeflux,
     856            &epsilon, &g, &H0)) {
    846857    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: flux_function_central could not parse input arguments");
    847858    return NULL;
    … …  
    850861 
    851862  _flux_function_central((double*) ql -> data,
    852                         (double*) qr -> data,
    853                         zl,
    854                          zr,                                           
    855                         ((double*) normal -> data)[0],
    856                          ((double*) normal -> data)[1],                 
    857                         epsilon, H0, g,
    858                         (double*) edgeflux -> data,
    859                         &max_speed);
     863            (double*) qr -> data,
     864            zl,
     865             zr,                       
     866            ((double*) normal -> data)[0],
     867             ((double*) normal -> data)[1],         
     868            epsilon, H0, g,
     869            (double*) edgeflux -> data,
     870            &max_speed);
    860871 
    861872  return Py_BuildValue("d", max_speed); 
    … …  
    877888
    878889  if (!PyArg_ParseTuple(args, "dOOOOO",
    879                         &g, &h, &v, &x,
    880                         &xmom, &ymom)) {
     890            &g, &h, &v, &x,
     891            &xmom, &ymom)) {
    881892    //&epsilon)) {
    882893    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: gravity could not parse input arguments");
    … …  
    936947
    937948  if (!PyArg_ParseTuple(args, "ddOOOOOOO",
    938                         &g, &eps, &w, &z, &uh, &vh, &eta,
    939                         &xmom, &ymom)) {
     949            &g, &eps, &w, &z, &uh, &vh, &eta,
     950            &xmom, &ymom)) {
    940951    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: manning_friction could not parse input arguments");
    941952    return NULL;
    … …  
    945956  N = w -> dimensions[0];
    946957  _manning_friction(g, eps, N,
    947                     (double*) w -> data,
    948                     (double*) z -> data,
    949                     (double*) uh -> data,
    950                     (double*) vh -> data,
    951                     (double*) eta -> data,
    952                     (double*) xmom -> data,
    953                     (double*) ymom -> data);
     958            (double*) w -> data,
     959            (double*) z -> data,
     960            (double*) uh -> data,
     961            (double*) vh -> data,
     962            (double*) eta -> data,
     963            (double*) xmom -> data,
     964            (double*) ymom -> data);
    954965
    955966  return Py_BuildValue("");
    … …  
    969980
    970981  if (!PyArg_ParseTuple(args, "ddOOOOOOO",
    971                         &g, &eps, &w, &z, &uh, &vh, &eta,
    972                         &xmom, &ymom))
     982            &g, &eps, &w, &z, &uh, &vh, &eta,
     983            &xmom, &ymom))
    973984    return NULL;
    974985
    975986  N = w -> dimensions[0];
    976987  _manning_friction_explicit(g, eps, N,
    977                     (double*) w -> data,
    978                     (double*) z -> data,
    979                     (double*) uh -> data,
    980                     (double*) vh -> data,
    981                     (double*) eta -> data,
    982                     (double*) xmom -> data,
    983                     (double*) ymom -> data);
     988            (double*) w -> data,
     989            (double*) z -> data,
     990            (double*) uh -> data,
     991            (double*) vh -> data,
     992            (double*) eta -> data,
     993            (double*) xmom -> data,
     994            (double*) ymom -> data);
    984995
    985996  return Py_BuildValue("");
    … …  
    9911002// Computational routine
    9921003int _extrapolate_second_order_sw(int number_of_elements,
    993                                   double epsilon,
    994                                   double minimum_allowed_height,
    995                                   double beta_w,
    996                                   double beta_w_dry,
    997                                   double beta_uh,
    998                                   double beta_uh_dry,
    999                                   double beta_vh,
    1000                                   double beta_vh_dry,
    1001                                   long* surrogate_neighbours,
    1002                                   long* number_of_boundaries,
    1003                                   double* centroid_coordinates,
    1004                                   double* stage_centroid_values,
    1005                                   double* xmom_centroid_values,
    1006                                   double* ymom_centroid_values,
    1007                                   double* elevation_centroid_values,
    1008                                   double* vertex_coordinates,
    1009                                   double* stage_vertex_values,
    1010                                   double* xmom_vertex_values,
    1011                                   double* ymom_vertex_values,
    1012                                   double* elevation_vertex_values,
    1013                                   int optimise_dry_cells) {
    1014                                  
    1015                                  
     1004                  double epsilon,
     1005                  double minimum_allowed_height,
     1006                  double beta_w,
     1007                  double beta_w_dry,
     1008                  double beta_uh,
     1009                  double beta_uh_dry,
     1010                  double beta_vh,
     1011                  double beta_vh_dry,
     1012                  long* surrogate_neighbours,
     1013                  long* number_of_boundaries,
     1014                  double* centroid_coordinates,
     1015                  double* stage_centroid_values,
     1016                  double* xmom_centroid_values,
     1017                  double* ymom_centroid_values,
     1018                  double* elevation_centroid_values,
     1019                  double* vertex_coordinates,
     1020                  double* stage_vertex_values,
     1021                  double* xmom_vertex_values,
     1022                  double* ymom_vertex_values,
     1023                  double* elevation_vertex_values,
     1024                  int optimise_dry_cells) {
     1025                 
     1026                 
    10161027
    10171028  // Local variables
    … …  
    10231034  double hc, h0, h1, h2, beta_tmp, hfactor;
    10241035
    1025        
     1036   
    10261037  for (k=0; k<number_of_elements; k++) {
    10271038    k3=k*3;
    … …  
    11081119      // the triangle might be flat or clockwise:
    11091120      if (area2<=0) {
    1110         PyErr_SetString(PyExc_RuntimeError,
    1111                         "shallow_water_ext.c: negative triangle area encountered");
    1112         return -1;
     1121    PyErr_SetString(PyExc_RuntimeError,
     1122            "shallow_water_ext.c: negative triangle area encountered");
     1123    return -1;
    11131124      } 
    11141125     
    … …  
    11231134      hfactor = 0.0;
    11241135      if (hmin > 0.001 ) {
    1125           hfactor = (hmin-0.001)/(hmin+0.004);
     1136      hfactor = (hmin-0.001)/(hmin+0.004);
    11261137      }
    11271138     
    11281139      if (optimise_dry_cells) {     
    1129         // Check if linear reconstruction is necessary for triangle k
    1130         // This check will exclude dry cells.
    1131 
    1132         hmax = max(h0,max(h1,h2));     
    1133         if (hmax < epsilon) {
    1134           continue;
    1135         }
     1140    // Check if linear reconstruction is necessary for triangle k
     1141    // This check will exclude dry cells.
     1142
     1143    hmax = max(h0,max(h1,h2));     
     1144    if (hmax < epsilon) {
     1145      continue;
     1146    }
    11361147      }
    11371148
    … …  
    11721183      //if (hmin>minimum_allowed_height)
    11731184      beta_tmp = beta_w_dry + (beta_w - beta_w_dry) * hfactor;
    1174        
     1185   
    11751186      //printf("min_alled_height = %f\n",minimum_allowed_height);
    11761187      //printf("hmin = %f\n",hmin);
    … …  
    11811192     
    11821193      for (i=0;i<3;i++)
    1183         stage_vertex_values[k3+i]=stage_centroid_values[k]+dqv[i];
     1194    stage_vertex_values[k3+i]=stage_centroid_values[k]+dqv[i];
    11841195     
    11851196     
    … …  
    12221233
    12231234      for (i=0;i<3;i++)
    1224         xmom_vertex_values[k3+i]=xmom_centroid_values[k]+dqv[i];
     1235    xmom_vertex_values[k3+i]=xmom_centroid_values[k]+dqv[i];
    12251236     
    12261237     
    … …  
    12591270      //if (hmin<minimum_allowed_height)
    12601271      //beta_tmp = beta_vh_dry;
    1261       beta_tmp = beta_vh_dry + (beta_vh - beta_vh_dry) * hfactor;       
     1272      beta_tmp = beta_vh_dry + (beta_vh - beta_vh_dry) * hfactor;   
    12621273
    12631274      // Limit the gradient
    … …  
    12651276     
    12661277      for (i=0;i<3;i++)
    1267         ymom_vertex_values[k3+i]=ymom_centroid_values[k]+dqv[i];
    1268        
     1278    ymom_vertex_values[k3+i]=ymom_centroid_values[k]+dqv[i];
     1279   
    12691280    } // End number_of_boundaries <=1
    12701281    else{
    … …  
    12781289      // Find the only internal neighbour (k1?)
    12791290      for (k2=k3;k2<k3+3;k2++){
    1280         // Find internal neighbour of triangle k     
    1281         // k2 indexes the edges of triangle k   
    1282      
    1283         if (surrogate_neighbours[k2]!=k)
    1284           break;
     1291    // Find internal neighbour of triangle k     
     1292    // k2 indexes the edges of triangle k   
     1293     
     1294    if (surrogate_neighbours[k2]!=k)
     1295      break;
    12851296      }
    12861297     
    12871298      if ((k2==k3+3)) {
    1288         // If we didn't find an internal neighbour
    1289         PyErr_SetString(PyExc_RuntimeError,
    1290                         "shallow_water_ext.c: Internal neighbour not found");     
    1291         return -1;
     1299    // If we didn't find an internal neighbour
     1300    PyErr_SetString(PyExc_RuntimeError,
     1301            "shallow_water_ext.c: Internal neighbour not found");     
     1302    return -1;
    12921303      }
    12931304     
    … …  
    13351346      // Now limit the jumps
    13361347      if (dq1>=0.0){
    1337         qmin=0.0;
    1338         qmax=dq1;
     1348    qmin=0.0;
     1349    qmax=dq1;
    13391350      }
    13401351      else{
    1341         qmin=dq1;
    1342         qmax=0.0;
     1352    qmin=dq1;
     1353    qmax=0.0;
    13431354      }
    13441355     
    … …  
    13471358     
    13481359      for (i=0;i<3;i++)
    1349         stage_vertex_values[k3+i]=stage_centroid_values[k]+dqv[i];
     1360    stage_vertex_values[k3+i]=stage_centroid_values[k]+dqv[i];
    13501361     
    13511362     
    … …  
    13691380      // Now limit the jumps
    13701381      if (dq1>=0.0){
    1371         qmin=0.0;
    1372         qmax=dq1;
     1382    qmin=0.0;
     1383    qmax=dq1;
    13731384      }
    13741385      else{
    1375         qmin=dq1;
    1376         qmax=0.0;
     1386    qmin=dq1;
     1387    qmax=0.0;
    13771388      }
    13781389     
    … …  
    13811392     
    13821393      for (i=0;i<3;i++)
    1383         xmom_vertex_values[k3+i]=xmom_centroid_values[k]+dqv[i];
     1394    xmom_vertex_values[k3+i]=xmom_centroid_values[k]+dqv[i];
    13841395     
    13851396     
    … …  
    14031414      // Now limit the jumps
    14041415      if (dq1>=0.0){
    1405         qmin=0.0;
    1406         qmax=dq1;
     1416    qmin=0.0;
     1417    qmax=dq1;
    14071418      }
    14081419      else{
    1409         qmin=dq1;
    1410         qmax=0.0;
     1420    qmin=dq1;
     1421    qmax=0.0;
    14111422      }
    14121423     
    … …  
    14151426     
    14161427      for (i=0;i<3;i++)
    1417         ymom_vertex_values[k3+i]=ymom_centroid_values[k]+dqv[i];
    1418        
     1428    ymom_vertex_values[k3+i]=ymom_centroid_values[k]+dqv[i];
     1429   
    14191430    } // else [number_of_boundaries==2]
    14201431  } // for k=0 to number_of_elements-1
    14211432 
    14221433  return 0;
    1423 }                       
     1434}           
    14241435
    14251436
    … …  
    14561467    Post conditions:
    14571468            The vertices of each triangle have values from a
    1458             limited linear reconstruction
    1459             based on centroid values
     1469        limited linear reconstruction
     1470        based on centroid values
    14601471
    14611472  */
    … …  
    14841495  // Convert Python arguments to C
    14851496  if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOi",
    1486                         &domain,
    1487                         &surrogate_neighbours,
    1488                         &number_of_boundaries,
    1489                         &centroid_coordinates,
    1490                         &stage_centroid_values,
    1491                         &xmom_centroid_values,
    1492                         &ymom_centroid_values,
    1493                         &elevation_centroid_values,
    1494                         &vertex_coordinates,
    1495                         &stage_vertex_values,
    1496                         &xmom_vertex_values,
    1497                         &ymom_vertex_values,
    1498                         &elevation_vertex_values,
    1499                         &optimise_dry_cells)) {                 
    1500                        
     1497            &domain,
     1498            &surrogate_neighbours,
     1499            &number_of_boundaries,
     1500            &centroid_coordinates,
     1501            &stage_centroid_values,
     1502            &xmom_centroid_values,
     1503            &ymom_centroid_values,
     1504            &elevation_centroid_values,
     1505            &vertex_coordinates,
     1506            &stage_vertex_values,
     1507            &xmom_vertex_values,
     1508            &ymom_vertex_values,
     1509            &elevation_vertex_values,
     1510            &optimise_dry_cells)) {         
     1511           
    15011512    PyErr_SetString(PyExc_RuntimeError,
    1502                     "Input arguments to extrapolate_second_order_sw failed");
     1513            "Input arguments to extrapolate_second_order_sw failed");
    15031514    return NULL;
    15041515  }
    … …  
    15221533  // Call underlying computational routine
    15231534  e = _extrapolate_second_order_sw(number_of_elements,
    1524                                    epsilon,
    1525                                    minimum_allowed_height,
    1526                                    beta_w,
    1527                                    beta_w_dry,
    1528                                    beta_uh,
    1529                                    beta_uh_dry,
    1530                                    beta_vh,
    1531                                    beta_vh_dry,
    1532                                    (long*) surrogate_neighbours -> data,
    1533                                    (long*) number_of_boundaries -> data,
    1534                                    (double*) centroid_coordinates -> data,
    1535                                    (double*) stage_centroid_values -> data,
    1536                                    (double*) xmom_centroid_values -> data,
    1537                                    (double*) ymom_centroid_values -> data,
    1538                                    (double*) elevation_centroid_values -> data,
    1539                                    (double*) vertex_coordinates -> data,
    1540                                    (double*) stage_vertex_values -> data,
    1541                                    (double*) xmom_vertex_values -> data,
    1542                                    (double*) ymom_vertex_values -> data,
    1543                                    (double*) elevation_vertex_values -> data,
    1544                                    optimise_dry_cells);
     1535                   epsilon,
     1536                   minimum_allowed_height,
     1537                   beta_w,
     1538                   beta_w_dry,
     1539                   beta_uh,
     1540                   beta_uh_dry,
     1541                   beta_vh,
     1542                   beta_vh_dry,
     1543                   (long*) surrogate_neighbours -> data,
     1544                   (long*) number_of_boundaries -> data,
     1545                   (double*) centroid_coordinates -> data,
     1546                   (double*) stage_centroid_values -> data,
     1547                   (double*) xmom_centroid_values -> data,
     1548                   (double*) ymom_centroid_values -> data,
     1549                   (double*) elevation_centroid_values -> data,
     1550                   (double*) vertex_coordinates -> data,
     1551                   (double*) stage_vertex_values -> data,
     1552                   (double*) xmom_vertex_values -> data,
     1553                   (double*) ymom_vertex_values -> data,
     1554                   (double*) elevation_vertex_values -> data,
     1555                   optimise_dry_cells);
    15451556  if (e == -1) {
    15461557    // Use error string set inside computational routine
    15471558    return NULL;
    1548   }                              
     1559  }                  
    15491560 
    15501561 
    … …  
    15741585  // Convert Python arguments to C
    15751586  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OO|i", argnames,
    1576                                    &Q, &Normal, &direction)) {
     1587                   &Q, &Normal, &direction)) {
    15771588    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: rotate could not parse input arguments");
    15781589    return NULL;
    … …  
    16191630// Computational function for flux computation
    16201631double _compute_fluxes_central(int number_of_elements,
    1621                                double timestep,
    1622                                double epsilon,
    1623                                double H0,
    1624                                double g,
    1625                                long* neighbours,
    1626                                long* neighbour_edges,
    1627                                double* normals,
    1628                                double* edgelengths,
    1629                                double* radii,
    1630                                double* areas,
    1631                                long* tri_full_flag,
    1632                                double* stage_edge_values,
    1633                                double* xmom_edge_values,
    1634                                double* ymom_edge_values,
    1635                                double* bed_edge_values,
    1636                                double* stage_boundary_values,
    1637                                double* xmom_boundary_values,
    1638                                double* ymom_boundary_values,
    1639                                double* stage_explicit_update,
    1640                                double* xmom_explicit_update,
    1641                                double* ymom_explicit_update,
    1642                                long* already_computed_flux,
    1643                                double* max_speed_array,
    1644                                int optimise_dry_cells) {
    1645                                
     1632                               double timestep,
     1633                               double epsilon,
     1634                               double H0,
     1635                               double g,
     1636                               long* neighbours,
     1637                               long* neighbour_edges,
     1638                               double* normals,
     1639                               double* edgelengths,
     1640                               double* radii,
     1641                               double* areas,
     1642                               long* tri_full_flag,
     1643                               double* stage_edge_values,
     1644                               double* xmom_edge_values,
     1645                               double* ymom_edge_values,
     1646                               double* bed_edge_values,
     1647                               double* stage_boundary_values,
     1648                               double* xmom_boundary_values,
     1649                               double* ymom_boundary_values,
     1650                               double* stage_explicit_update,
     1651                               double* xmom_explicit_update,
     1652                               double* ymom_explicit_update,
     1653                               long* already_computed_flux,
     1654                               double* max_speed_array,
     1655                               int optimise_dry_cells)
     1656{                   
    16461657  // Local variables
    1647   double max_speed, length, area, zl, zr;
     1658  double max_speed, length, inv_area, zl, zr;
    16481659  int k, i, m, n;
    16491660  int ki, nm=0, ki2; // Index shorthands
    … …  
    16521663  double ql[3], qr[3], edgeflux[3]; // Work array for summing up fluxes
    16531664
    1654   static long call=1; // Static local variable flagging already computed flux
    1655                                
    1656        
     1665  static long call = 1; // Static local variable flagging already computed flux
     1666                   
    16571667  // Start computation
    16581668  call++; // Flag 'id' of flux calculation for this timestep
    16591669
    1660    
    16611670  // Set explicit_update to zero for all conserved_quantities.
    16621671  // This assumes compute_fluxes called before forcing terms
    … …  
    16651674  memset((char*) ymom_explicit_update, 0, number_of_elements*sizeof(double));   
    16661675 
    1667 
    16681676  // For all triangles
    1669   for (k=0; k<number_of_elements; k++) {
    1670    
     1677  for (k = 0; k < number_of_elements; k++)
     1678  { 
    16711679    // Loop through neighbours and compute edge flux for each 
    1672     for (i=0; i<3; i++) {
    1673       ki = k*3+i; // Linear index to edge i of triangle k
     1680    for (i = 0; i < 3; i++)
     1681    {
     1682      ki = k*3 + i; // Linear index to edge i of triangle k
    16741683     
    16751684      if (already_computed_flux[ki] == call)
     1685      {
    16761686        // We've already computed the flux across this edge
    1677         continue;
    1678        
     1687        continue;
     1688      }
     1689   
    16791690      // Get left hand side values from triangle k, edge i
    16801691      ql[0] = stage_edge_values[ki];
    … …  
    16861697      // or from boundary array (Quantities at neighbour on nearest face).
    16871698      n = neighbours[ki];
    1688       if (n < 0) {
     1699      if (n < 0)
     1700      {
    16891701        // Neighbour is a boundary condition
    1690         m = -n-1; // Convert negative flag to boundary index
    1691        
    1692         qr[0] = stage_boundary_values[m];
    1693         qr[1] = xmom_boundary_values[m];
    1694         qr[2] = ymom_boundary_values[m];
    1695         zr = zl; // Extend bed elevation to boundary
    1696       } else {
     1702        m = -n - 1; // Convert negative flag to boundary index
     1703   
     1704        qr[0] = stage_boundary_values[m];
     1705        qr[1] = xmom_boundary_values[m];
     1706        qr[2] = ymom_boundary_values[m];
     1707        zr = zl; // Extend bed elevation to boundary
     1708      }
     1709      else
     1710      {
    16971711        // Neighbour is a real triangle
    1698         m = neighbour_edges[ki];
    1699         nm = n*3+m; // Linear index (triangle n, edge m)
    1700        
    1701         qr[0] = stage_edge_values[nm];
    1702         qr[1] = xmom_edge_values[nm];
    1703         qr[2] = ymom_edge_values[nm];
    1704         zr = bed_edge_values[nm];
     1712        m = neighbour_edges[ki];
     1713        nm = n*3 + m; // Linear index (triangle n, edge m)
     1714       
     1715        qr[0] = stage_edge_values[nm];
     1716        qr[1] = xmom_edge_values[nm];
     1717        qr[2] = ymom_edge_values[nm];
     1718        zr = bed_edge_values[nm];
    17051719      }
    17061720     
    17071721      // Now we have values for this edge - both from left and right side.
    17081722     
    1709       if (optimise_dry_cells) {     
    1710         // Check if flux calculation is necessary across this edge
    1711         // This check will exclude dry cells.
    1712         // This will also optimise cases where zl != zr as
    1713         // long as both are dry
    1714 
    1715         if ( fabs(ql[0] - zl) < epsilon &&
    1716              fabs(qr[0] - zr) < epsilon ) {
    1717           // Cell boundary is dry
    1718          
    1719           already_computed_flux[ki] = call; // #k Done 
    1720           if (n>=0)
    1721             already_computed_flux[nm] = call; // #n Done
    1722        
    1723           max_speed = 0.0;
    1724           continue;
    1725         }
    1726       }
    1727    
     1723      if (optimise_dry_cells)
     1724      {     
     1725        // Check if flux calculation is necessary across this edge
     1726        // This check will exclude dry cells.
     1727        // This will also optimise cases where zl != zr as
     1728        // long as both are dry
     1729
     1730        if (fabs(ql[0] - zl) < epsilon &&
     1731            fabs(qr[0] - zr) < epsilon)
     1732        {
     1733          // Cell boundary is dry
     1734         
     1735          already_computed_flux[ki] = call; // #k Done 
     1736          if (n >= 0)
     1737          {
     1738            already_computed_flux[nm] = call; // #n Done
     1739          }
     1740       
     1741          max_speed = 0.0;
     1742          continue;
     1743        }
     1744      }
    17281745     
    17291746      // Outward pointing normal vector (domain.normals[k, 2*i:2*i+2])
    … …  
    17321749      // Edge flux computation (triangle k, edge i)
    17331750      _flux_function_central(ql, qr, zl, zr,
    1734                              normals[ki2], normals[ki2+1],
    1735                              epsilon, H0, g,
    1736                              edgeflux, &max_speed);
     1751                             normals[ki2], normals[ki2+1],
     1752                             epsilon, H0, g,
     1753                             edgeflux, &max_speed);
    17371754     
    17381755     
    … …  
    17531770     
    17541771      // Update neighbour n with same flux but reversed sign
    1755       if (n>=0) {
    1756         stage_explicit_update[n] += edgeflux[0];
    1757         xmom_explicit_update[n] += edgeflux[1];
    1758         ymom_explicit_update[n] += edgeflux[2];
    1759        
    1760         already_computed_flux[nm] = call; // #n Done
    1761       }
    1762 
     1772      if (n >= 0)
     1773      {
     1774        stage_explicit_update[n] += edgeflux[0];
     1775        xmom_explicit_update[n] += edgeflux[1];
     1776        ymom_explicit_update[n] += edgeflux[2];
     1777       
     1778        already_computed_flux[nm] = call; // #n Done
     1779      }
    17631780
    17641781      // Update timestep based on edge i and possibly neighbour n
    1765       if (tri_full_flag[k] == 1) {
    1766         if (max_speed > epsilon) {
    1767        
    1768           // Apply CFL condition for triangles joining this edge (triangle k and triangle n)
    1769          
    1770           // CFL for triangle k
    1771           timestep = min(timestep, radii[k]/max_speed);
    1772          
    1773           if (n>=0)
    1774             // Apply CFL condition for neigbour n (which is on the ith edge of triangle k)
    1775             timestep = min(timestep, radii[n]/max_speed);
    1776          
    1777           // Ted Rigby's suggested less conservative version
    1778           //if (n>=0) {             
    1779           //  timestep = min(timestep, (radii[k]+radii[n])/max_speed);
    1780           //} else {
    1781           //  timestep = min(timestep, radii[k]/max_speed);
    1782           // }
    1783         }
     1782      if (tri_full_flag[k] == 1)
     1783      {
     1784        if (max_speed > epsilon)
     1785        {
     1786          // Apply CFL condition for triangles joining this edge (triangle k and triangle n)
     1787         
     1788          // CFL for triangle k
     1789          timestep = min(timestep, radii[k]/max_speed);
     1790         
     1791          if (n >= 0)
     1792          {
     1793            // Apply CFL condition for neigbour n (which is on the ith edge of triangle k)
     1794            timestep = min(timestep, radii[n]/max_speed);
     1795          }
     1796
     1797          // Ted Rigby's suggested less conservative version
     1798          //if (n>=0) {         
     1799          //  timestep = min(timestep, (radii[k]+radii[n])/max_speed);
     1800          //} else {
     1801          //  timestep = min(timestep, radii[k]/max_speed);
     1802          // }
     1803        }
    17841804      }
    17851805     
    … …  
    17891809    // Normalise triangle k by area and store for when all conserved
    17901810    // quantities get updated
    1791     area = areas[k];
    1792     stage_explicit_update[k] /= area;
    1793     xmom_explicit_update[k] /= area;
    1794     ymom_explicit_update[k] /= area;
     1811    inv_area = 1.0/areas[k];
     1812    stage_explicit_update[k] *= inv_area;
     1813    xmom_explicit_update[k] *= inv_area;
     1814    ymom_explicit_update[k] *= inv_area;
    17951815   
    17961816   
    … …  
    17991819   
    18001820  } // End triangle k
    1801        
    1802                                
    1803                                
     1821             
    18041822  return timestep;
    18051823}
    … …  
    18361854
    18371855    PyObject
    1838         *domain,
    1839         *stage,
    1840         *xmom,
    1841         *ymom,
    1842         *bed;
     1856    *domain,
     1857    *stage,
     1858    *xmom,
     1859    *ymom,
     1860    *bed;
    18431861
    18441862    PyArrayObject
    1845         *neighbours,
    1846         *neighbour_edges,
    1847         *normals,
    1848         *edgelengths,
    1849         *radii,
    1850         *areas,
    1851         *tri_full_flag,
    1852         *stage_edge_values,
    1853         *xmom_edge_values,
    1854         *ymom_edge_values,
    1855         *bed_edge_values,
    1856         *stage_boundary_values,
    1857         *xmom_boundary_values,
    1858         *ymom_boundary_values,
    1859         *stage_explicit_update,
    1860         *xmom_explicit_update,
    1861         *ymom_explicit_update,
    1862         *already_computed_flux, //Tracks whether the flux across an edge has already been computed
    1863         *max_speed_array; //Keeps track of max speeds for each triangle
     1863    *neighbours,
     1864    *neighbour_edges,
     1865    *normals,
     1866    *edgelengths,
     1867    *radii,
     1868    *areas,
     1869    *tri_full_flag,
     1870    *stage_edge_values,
     1871    *xmom_edge_values,
     1872    *ymom_edge_values,
     1873    *bed_edge_values,
     1874    *stage_boundary_values,
     1875    *xmom_boundary_values,
     1876    *ymom_boundary_values,
     1877    *stage_explicit_update,
     1878    *xmom_explicit_update,
     1879    *ymom_explicit_update,
     1880    *already_computed_flux, //Tracks whether the flux across an edge has already been computed
     1881    *max_speed_array; //Keeps track of max speeds for each triangle
    18641882
    18651883   
    … …  
    18691887    // Convert Python arguments to C
    18701888    if (!PyArg_ParseTuple(args, "dOOOO", &timestep, &domain, &stage, &xmom, &ymom, &bed )) {
    1871         PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
    1872         return NULL;
     1889    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
     1890    return NULL;
    18731891    }
    18741892
    … …  
    19071925  // the explicit update arrays
    19081926  timestep = _compute_fluxes_central(number_of_elements,
    1909                                      timestep,
    1910                                      epsilon,
    1911                                      H0,
    1912                                      g,
    1913                                      (long*) neighbours -> data,
    1914                                      (long*) neighbour_edges -> data,
    1915                                      (double*) normals -> data,
    1916                                      (double*) edgelengths -> data,
    1917                                      (double*) radii -> data,
    1918                                      (double*) areas -> data,
    1919                                      (long*) tri_full_flag -> data,
    1920                                      (double*) stage_edge_values -> data,
    1921                                      (double*) xmom_edge_values -> data,
    1922                                      (double*) ymom_edge_values -> data,
    1923                                      (double*) bed_edge_values -> data,
    1924                                      (double*) stage_boundary_values -> data,
    1925                                      (double*) xmom_boundary_values -> data,
    1926                                      (double*) ymom_boundary_values -> data,
    1927                                      (double*) stage_explicit_update -> data,
    1928                                      (double*) xmom_explicit_update -> data,
    1929                                      (double*) ymom_explicit_update -> data,
    1930                                      (long*) already_computed_flux -> data,
    1931                                      (double*) max_speed_array -> data,
    1932                                      optimise_dry_cells);
     1927                     timestep,
     1928                     epsilon,
     1929                     H0,
     1930                     g,
     1931                     (long*) neighbours -> data,
     1932                     (long*) neighbour_edges -> data,
     1933                     (double*) normals -> data,
     1934                     (double*) edgelengths -> data,
     1935                     (double*) radii -> data,
     1936                     (double*) areas -> data,
     1937                     (long*) tri_full_flag -> data,
     1938                     (double*) stage_edge_values -> data,
     1939                     (double*) xmom_edge_values -> data,
     1940                     (double*) ymom_edge_values -> data,
     1941                     (double*) bed_edge_values -> data,
     1942                     (double*) stage_boundary_values -> data,
     1943                     (double*) xmom_boundary_values -> data,
     1944                     (double*) ymom_boundary_values -> data,
     1945                     (double*) stage_explicit_update -> data,
     1946                     (double*) xmom_explicit_update -> data,
     1947                     (double*) ymom_explicit_update -> data,
     1948                     (long*) already_computed_flux -> data,
     1949                     (double*) max_speed_array -> data,
     1950                     optimise_dry_cells);
    19331951
    19341952  Py_DECREF(neighbours);
    … …  
    19801998    domain.timestep = compute_fluxes(timestep,
    19811999                                     domain.epsilon,
    1982                                      domain.H0,
     2000                     domain.H0,
    19832001                                     domain.g,
    19842002                                     domain.neighbours,
    … …  
    20002018                                     Ymom.explicit_update,
    20012019                                     already_computed_flux,
    2002                                      optimise_dry_cells)                                     
     2020                     optimise_dry_cells)                     
    20032021
    20042022
    … …  
    20332051  // Convert Python arguments to C
    20342052  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOi",
    2035                         &timestep,
    2036                         &epsilon,
    2037                         &H0,
    2038                         &g,
    2039                         &neighbours,
    2040                         &neighbour_edges,
    2041                         &normals,
    2042                         &edgelengths, &radii, &areas,
    2043                         &tri_full_flag,
    2044                         &stage_edge_values,
    2045                         &xmom_edge_values,
    2046                         &ymom_edge_values,
    2047                         &bed_edge_values,
    2048                         &stage_boundary_values,
    2049                         &xmom_boundary_values,
    2050                         &ymom_boundary_values,
    2051                         &stage_explicit_update,
    2052                         &xmom_explicit_update,
    2053                         &ymom_explicit_update,
    2054                         &already_computed_flux,
    2055                         &max_speed_array,
    2056                         &optimise_dry_cells)) {
     2053            &timestep,
     2054            &epsilon,
     2055            &H0,
     2056            &g,
     2057            &neighbours,
     2058            &neighbour_edges,
     2059            &normals,
     2060            &edgelengths, &radii, &areas,
     2061            &tri_full_flag,
     2062            &stage_edge_values,
     2063            &xmom_edge_values,
     2064            &ymom_edge_values,
     2065            &bed_edge_values,
     2066            &stage_boundary_values,
     2067            &xmom_boundary_values,
     2068            &ymom_boundary_values,
     2069            &stage_explicit_update,
     2070            &xmom_explicit_update,
     2071            &ymom_explicit_update,
     2072            &already_computed_flux,
     2073            &max_speed_array,
     2074            &optimise_dry_cells)) {
    20572075    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
    20582076    return NULL;
    … …  
    20652083  // the explicit update arrays
    20662084  timestep = _compute_fluxes_central(number_of_elements,
    2067                                      timestep,
    2068                                      epsilon,
    2069                                      H0,
    2070                                      g,
    2071                                      (long*) neighbours -> data,
    2072                                      (long*) neighbour_edges -> data,
    2073                                      (double*) normals -> data,
    2074                                      (double*) edgelengths -> data,
    2075                                      (double*) radii -> data,
    2076                                      (double*) areas -> data,
    2077                                      (long*) tri_full_flag -> data,
    2078                                      (double*) stage_edge_values -> data,
    2079                                      (double*) xmom_edge_values -> data,
    2080                                      (double*) ymom_edge_values -> data,
    2081                                      (double*) bed_edge_values -> data,
    2082                                      (double*) stage_boundary_values -> data,
    2083                                      (double*) xmom_boundary_values -> data,
    2084                                      (double*) ymom_boundary_values -> data,
    2085                                      (double*) stage_explicit_update -> data,
    2086                                      (double*) xmom_explicit_update -> data,
    2087                                      (double*) ymom_explicit_update -> data,
    2088                                      (long*) already_computed_flux -> data,
    2089                                      (double*) max_speed_array -> data,
    2090                                      optimise_dry_cells);
     2085                     timestep,
     2086                     epsilon,
     2087                     H0,
     2088                     g,
     2089                     (long*) neighbours -> data,
     2090                     (long*) neighbour_edges -> data,
     2091                     (double*) normals -> data,
     2092                     (double*) edgelengths -> data,
     2093                     (double*) radii -> data,
     2094                     (double*) areas -> data,
     2095                     (long*) tri_full_flag -> data,
     2096                     (double*) stage_edge_values -> data,
     2097                     (double*) xmom_edge_values -> data,
     2098                     (double*) ymom_edge_values -> data,
     2099                     (double*) bed_edge_values -> data,
     2100                     (double*) stage_boundary_values -> data,
     2101                     (double*) xmom_boundary_values -> data,
     2102                     (double*) ymom_boundary_values -> data,
     2103                     (double*) stage_explicit_update -> data,
     2104                     (double*) xmom_explicit_update -> data,
     2105                     (double*) ymom_explicit_update -> data,
     2106                     (long*) already_computed_flux -> data,
     2107                     (double*) max_speed_array -> data,
     2108                     optimise_dry_cells);
    20912109 
    20922110  // Return updated flux timestep
    … …  
    21752193  // Convert Python arguments to C
    21762194  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOO",
    2177                         &timestep,
    2178                         &epsilon,
    2179                         &H0,
    2180                         &g,
    2181                         &neighbours,
    2182                         &neighbour_edges,
    2183                         &normals,
    2184                         &edgelengths, &radii, &areas,
    2185                         &tri_full_flag,
    2186                         &stage_edge_values,
    2187                         &xmom_edge_values,
    2188                         &ymom_edge_values,
    2189                         &bed_edge_values,
    2190                         &stage_boundary_values,
    2191                         &xmom_boundary_values,
    2192                         &ymom_boundary_values,
    2193                         &stage_explicit_update,
    2194                         &xmom_explicit_update,
    2195                         &ymom_explicit_update,
    2196                         &already_computed_flux)) {
     2195            &timestep,
     2196            &epsilon,
     2197            &H0,
     2198            &g,
     2199            &neighbours,
     2200            &neighbour_edges,
     2201            &normals,
     2202            &edgelengths, &radii, &areas,
     2203            &tri_full_flag,
     2204            &stage_edge_values,
     2205            &xmom_edge_values,
     2206            &ymom_edge_values,
     2207            &bed_edge_values,
     2208            &stage_boundary_values,
     2209            &xmom_boundary_values,
     2210            &ymom_boundary_values,
     2211            &stage_explicit_update,
     2212            &xmom_explicit_update,
     2213            &ymom_explicit_update,
     2214            &already_computed_flux)) {
    21972215    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
    21982216    return NULL;
    … …  
    22122230      ki = k*3+i;
    22132231      if (((long *) already_computed_flux->data)[ki]==call)//we've already computed the flux across this edge
    2214         continue;
     2232    continue;
    22152233      ql[0] = ((double *) stage_edge_values -> data)[ki];
    22162234      ql[1] = ((double *) xmom_edge_values -> data)[ki];
    … …  
    22212239      n = ((long *) neighbours -> data)[ki];
    22222240      if (n < 0) {
    2223         m = -n-1; //Convert negative flag to index
    2224         qr[0] = ((double *) stage_boundary_values -> data)[m];
    2225         qr[1] = ((double *) xmom_boundary_values -> data)[m];
    2226         qr[2] = ((double *) ymom_boundary_values -> data)[m];
    2227         zr = zl; //Extend bed elevation to boundary
     2241    m = -n-1; //Convert negative flag to index
     2242    qr[0] = ((double *) stage_boundary_values -> data)[m];
     2243    qr[1] = ((double *) xmom_boundary_values -> data)[m];
     2244    qr[2] = ((double *) ymom_boundary_values -> data)[m];
     2245    zr = zl; //Extend bed elevation to boundary
    22282246      } else {
    2229         m = ((long *) neighbour_edges -> data)[ki];
    2230         nm = n*3+m;
    2231         qr[0] = ((double *) stage_edge_values -> data)[nm];
    2232         qr[1] = ((double *) xmom_edge_values -> data)[nm];
    2233         qr[2] = ((double *) ymom_edge_values -> data)[nm];
    2234         zr =    ((double *) bed_edge_values -> data)[nm];
     2247    m = ((long *) neighbour_edges -> data)[ki];
     2248    nm = n*3+m;
     2249    qr[0] = ((double *) stage_edge_values -> data)[nm];
     2250    qr[1] = ((double *) xmom_edge_values -> data)[nm];
     2251    qr[2] = ((double *) ymom_edge_values -> data)[nm];
     2252    zr =    ((double *) bed_edge_values -> data)[nm];
    22352253      }
    22362254      // Outward pointing normal vector
    … …  
    22412259      //Edge flux computation
    22422260      flux_function_kinetic(ql, qr, zl, zr,
    2243                     normal[0], normal[1],
    2244                     epsilon, H0, g,
    2245                     edgeflux, &max_speed);
     2261            normal[0], normal[1],
     2262            epsilon, H0, g,
     2263            edgeflux, &max_speed);
    22462264      //update triangle k
    22472265      ((long *) already_computed_flux->data)[ki]=call;
    … …  
    22512269      //update the neighbour n
    22522270      if (n>=0){
    2253         ((long *) already_computed_flux->data)[nm]=call;
    2254         ((double *) stage_explicit_update -> data)[n] += edgeflux[0]*((double *) edgelengths -> data)[nm];
    2255         ((double *) xmom_explicit_update -> data)[n] += edgeflux[1]*((double *) edgelengths -> data)[nm];
    2256         ((double *) ymom_explicit_update -> data)[n] += edgeflux[2]*((double *) edgelengths -> data)[nm];
     2271    ((long *) already_computed_flux->data)[nm]=call;
     2272    ((double *) stage_explicit_update -> data)[n] += edgeflux[0]*((double *) edgelengths -> data)[nm];
     2273    ((double *) xmom_explicit_update -> data)[n] += edgeflux[1]*((double *) edgelengths -> data)[nm];
     2274    ((double *) ymom_explicit_update -> data)[n] += edgeflux[2]*((double *) edgelengths -> data)[nm];
    22572275      }
    22582276      ///for (j=0; j<3; j++) {
    2259         ///flux[j] -= edgeflux[j]*((double *) edgelengths -> data)[ki];
    2260         ///}
    2261         //Update timestep
    2262         //timestep = min(timestep, domain.radii[k]/max_speed)
    2263         //FIXME: SR Add parameter for CFL condition
     2277    ///flux[j] -= edgeflux[j]*((double *) edgelengths -> data)[ki];
     2278    ///}
     2279    //Update timestep
     2280    //timestep = min(timestep, domain.radii[k]/max_speed)
     2281    //FIXME: SR Add parameter for CFL condition
    22642282    if ( ((long *) tri_full_flag -> data)[k] == 1) {
    2265             if (max_speed > epsilon) {
    2266                 timestep = min(timestep, ((double *) radii -> data)[k]/max_speed);
    2267                 //maxspeed in flux_function is calculated as max(|u+a|,|u-a|)
    2268                 if (n>=0)
    2269                     timestep = min(timestep, ((double *) radii -> data)[n]/max_speed);
    2270             }
     2283        if (max_speed > epsilon) {
     2284            timestep = min(timestep, ((double *) radii -> data)[k]/max_speed);
     2285            //maxspeed in flux_function is calculated as max(|u+a|,|u-a|)
     2286            if (n>=0)
     2287                timestep = min(timestep, ((double *) radii -> data)[n]/max_speed);
     2288        }
    22712289    }
    22722290    } // end for i
    … …  
    22972315  // Convert Python arguments to C
    22982316  if (!PyArg_ParseTuple(args, "dddOOOO",
    2299                         &minimum_allowed_height,
    2300                         &maximum_allowed_speed,
    2301                         &epsilon,
    2302                         &wc, &zc, &xmomc, &ymomc)) {
     2317            &minimum_allowed_height,
     2318            &maximum_allowed_speed,
     2319            &epsilon,
     2320            &wc, &zc, &xmomc, &ymomc)) {
    23032321    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: protect could not parse input arguments");
    23042322    return NULL;
    … …  
    23082326
    23092327  _protect(N,
    2310            minimum_allowed_height,
    2311            maximum_allowed_speed,
    2312            epsilon,
    2313            (double*) wc -> data,
    2314            (double*) zc -> data,
    2315            (double*) xmomc -> data,
    2316            (double*) ymomc -> data);
     2328       minimum_allowed_height,
     2329       maximum_allowed_speed,
     2330       epsilon,
     2331       (double*) wc -> data,
     2332       (double*) zc -> data,
     2333       (double*) xmomc -> data,
     2334       (double*) ymomc -> data);
    23172335
    23182336  return Py_BuildValue("");
    … …  
    23562374  // Convert Python arguments to C
    23572375  if (!PyArg_ParseTuple(args, "OOOOOOOOOO",
    2358                         &domain,
    2359                         &wc, &zc,
    2360                         &wv, &zv, &hvbar,
    2361                         &xmomc, &ymomc, &xmomv, &ymomv)) {
     2376            &domain,
     2377            &wc, &zc,
     2378            &wv, &zv, &hvbar,
     2379            &xmomc, &ymomc, &xmomv, &ymomv)) {
    23622380    PyErr_SetString(PyExc_RuntimeError,
    2363                     "shallow_water_ext.c: balance_deep_and_shallow could not parse input arguments");
     2381            "shallow_water_ext.c: balance_deep_and_shallow could not parse input arguments");
    23642382    return NULL;
    23652383  } 
    2366          
    2367          
     2384     
     2385     
    23682386  // FIXME (Ole): I tested this without GetAttrString and got time down
    23692387  // marginally from 4.0s to 3.8s. Consider passing everything in
    … …  
    24102428  N = wc -> dimensions[0];
    24112429  _balance_deep_and_shallow(N,
    2412                             (double*) wc -> data,
    2413                             (double*) zc -> data,
    2414                             (double*) wv -> data,
    2415                             (double*) zv -> data,
    2416                             (double*) hvbar -> data,
    2417                             (double*) xmomc -> data,
    2418                             (double*) ymomc -> data,
    2419                             (double*) xmomv -> data,
    2420                             (double*) ymomv -> data,
    2421                             H0,
    2422                             (int) tight_slope_limiters,
    2423                             (int) use_centroid_velocities,                         
    2424                             alpha_balance);
     2430                (double*) wc -> data,
     2431                (double*) zc -> data,
     2432                (double*) wv -> data,
     2433                (double*) zv -> data,
     2434                (double*) hvbar -> data,
     2435                (double*) xmomc -> data,
     2436                (double*) ymomc -> data,
     2437                (double*) xmomv -> data,
     2438                (double*) ymomv -> data,
     2439                H0,
     2440                (int) tight_slope_limiters,
     2441                (int) use_centroid_velocities,             
     2442                alpha_balance);
    24252443
    24262444
    … …  
    24502468  // Convert Python arguments to C
    24512469  if (!PyArg_ParseTuple(args, "OOOOd",
    2452                         &xmom_update,
    2453                         &ymom_update,
    2454                         &s_vec, &phi_vec,
    2455                         &cw)) {
     2470            &xmom_update,
     2471            &ymom_update,
     2472            &s_vec, &phi_vec,
     2473            &cw)) {
    24562474    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: assign_windfield_values could not parse input arguments");
    24572475    return NULL;
    24582476  } 
    2459                        
     2477           
    24602478
    24612479  N = xmom_update -> dimensions[0];
    24622480
    24632481  _assign_wind_field_values(N,
    2464            (double*) xmom_update -> data,
    2465            (double*) ymom_update -> data,
    2466            (double*) s_vec -> data,
    2467            (double*) phi_vec -> data,
    2468            cw);
     2482       (double*) xmom_update -> data,
     2483       (double*) ymom_update -> data,
     2484       (double*) s_vec -> data,
     2485       (double*) phi_vec -> data,
     2486       cw);
    24692487
    24702488  return Py_BuildValue("");
Note: See TracChangeset for help on using the changeset viewer.