source: trunk/anuga_core/source/anuga/shallow_water/shallow_water_ext.c @ 8017

Last change on this file since 8017 was 8017, checked in by steve, 14 years ago

Added code to pickup if elevation is discontinuous. compute_fluxes now produces an error in this case

File size: 87.3 KB
Line 
1// Python - C extension module for shallow_water.py
2//
3// To compile (Python2.3):
4//  gcc -c domain_ext.c -I/usr/include/python2.3 -o domain_ext.o -Wall -O
5//  gcc -shared domain_ext.o  -o domain_ext.so
6//
7// or use python compile.py
8//
9// See the module shallow_water_domain.py for more documentation on
10// how to use this module
11//
12//
13// Ole Nielsen, GA 2004
14
15
16#include "Python.h"
17#include "numpy/arrayobject.h"
18#include "math.h"
19#include <stdio.h>
20#include "numpy_shim.h"
21
22// Shared code snippets
23#include "util_ext.h"
24
25
26
27
28
29const double pi = 3.14159265358979;
30
31// Computational function for rotation
32// Tried to inline, but no speedup was achieved 27th May 2009 (Ole)
33// static inline int _rotate(double *q, double n1, double n2) {
34int _rotate(double *q, double n1, double n2) {
35  /*Rotate the momentum component q (q[1], q[2])
36    from x,y coordinates to coordinates based on normal vector (n1, n2).
37
38    Result is returned in array 3x1 r
39    To rotate in opposite direction, call rotate with (q, n1, -n2)
40
41    Contents of q are changed by this function */
42
43
44  double q1, q2;
45
46  // Shorthands
47  q1 = q[1];  // uh momentum
48  q2 = q[2];  // vh momentum
49
50  // Rotate
51  q[1] =  n1*q1 + n2*q2;
52  q[2] = -n2*q1 + n1*q2;
53
54  return 0;
55}
56
57int find_qmin_and_qmax(double dq0, double dq1, double dq2, 
58               double *qmin, double *qmax){
59  // Considering the centroid of an FV triangle and the vertices of its
60  // auxiliary triangle, find
61  // qmin=min(q)-qc and qmax=max(q)-qc,
62  // where min(q) and max(q) are respectively min and max over the
63  // four values (at the centroid of the FV triangle and the auxiliary
64  // triangle vertices),
65  // and qc is the centroid
66  // dq0=q(vertex0)-q(centroid of FV triangle)
67  // dq1=q(vertex1)-q(vertex0)
68  // dq2=q(vertex2)-q(vertex0)
69 
70  if (dq0>=0.0){
71    if (dq1>=dq2){
72      if (dq1>=0.0)
73        *qmax=dq0+dq1;
74      else
75        *qmax=dq0;
76     
77      *qmin=dq0+dq2;
78      if (*qmin>=0.0) *qmin = 0.0;
79    }
80    else{// dq1<dq2
81      if (dq2>0)
82        *qmax=dq0+dq2;
83      else
84        *qmax=dq0;
85     
86      *qmin=dq0+dq1;   
87      if (*qmin>=0.0) *qmin=0.0;
88    }
89  }
90  else{//dq0<0
91    if (dq1<=dq2){
92      if (dq1<0.0)
93        *qmin=dq0+dq1;
94      else
95        *qmin=dq0;
96     
97      *qmax=dq0+dq2;   
98      if (*qmax<=0.0) *qmax=0.0;
99    }
100    else{// dq1>dq2
101      if (dq2<0.0)
102        *qmin=dq0+dq2;
103      else
104        *qmin=dq0;
105     
106      *qmax=dq0+dq1;
107      if (*qmax<=0.0) *qmax=0.0;
108    }
109  }
110  return 0;
111}
112
113int limit_gradient(double *dqv, double qmin, double qmax, double beta_w){
114  // Given provisional jumps dqv from the FV triangle centroid to its
115  // vertices and jumps qmin (qmax) between the centroid of the FV
116  // triangle and the minimum (maximum) of the values at the centroid of
117  // the FV triangle and the auxiliary triangle vertices,
118  // calculate a multiplicative factor phi by which the provisional
119  // vertex jumps are to be limited
120 
121  int i;
122  double r=1000.0, r0=1.0, phi=1.0;
123  static double TINY = 1.0e-100; // to avoid machine accuracy problems.
124  // FIXME: Perhaps use the epsilon used elsewhere.
125 
126  // Any provisional jump with magnitude < TINY does not contribute to
127  // the limiting process.
128  for (i=0;i<3;i++){
129    if (dqv[i]<-TINY)
130      r0=qmin/dqv[i];
131     
132    if (dqv[i]>TINY)
133      r0=qmax/dqv[i];
134     
135    r=min(r0,r);
136  }
137 
138  phi=min(r*beta_w,1.0);
139  //for (i=0;i<3;i++)
140  dqv[0]=dqv[0]*phi;
141  dqv[1]=dqv[1]*phi;
142  dqv[2]=dqv[2]*phi;
143   
144  return 0;
145}
146
147
148double compute_froude_number(double uh,
149                 double h, 
150                 double g,
151                 double epsilon) {
152             
153  // Compute Froude number; v/sqrt(gh)
154  // FIXME (Ole): Not currently in use
155   
156  double froude_number;
157 
158  //Compute Froude number (stability diagnostics)
159  if (h > epsilon) { 
160    froude_number = uh/sqrt(g*h)/h;
161  } else {
162    froude_number = 0.0;
163    // FIXME (Ole): What should it be when dry??
164  }
165 
166  return froude_number;
167}
168
169
170
171// Function to obtain speed from momentum and depth.
172// This is used by flux functions
173// Input parameters uh and h may be modified by this function.
174// Tried to inline, but no speedup was achieved 27th May 2009 (Ole)
175//static inline double _compute_speed(double *uh,
176double _compute_speed(double *uh, 
177                      double *h, 
178                      double epsilon, 
179                      double h0,
180                      double limiting_threshold) {
181 
182  double u;
183
184  if (*h < limiting_threshold) {   
185    // Apply limiting of speeds according to the ANUGA manual
186    if (*h < epsilon) {
187      *h = 0.0;  // Could have been negative
188      u = 0.0;
189    } else {
190      u = *uh/(*h + h0/ *h);   
191    }
192 
193
194    // Adjust momentum to be consistent with speed
195    *uh = u * *h;
196  } else {
197    // We are in deep water - no need for limiting
198    u = *uh/ *h;
199  }
200 
201  return u;
202}
203
204
205// Optimised squareroot computation
206//
207//See
208//http://www.lomont.org/Math/Papers/2003/InvSqrt.pdf
209//and http://mail.opensolaris.org/pipermail/tools-gcc/2005-August/000047.html
210float fast_squareroot_approximation(float number) {
211  float x; 
212  const float f = 1.5;
213
214  // Allow i and y to occupy the same memory
215  union
216  {
217    long i;
218    float y;
219  } u;
220 
221  // Good initial guess 
222  x = number * 0.5; 
223  u.= number;
224  u.= 0x5f3759df - ( u.i >> 1 );
225 
226  // Take a few iterations
227  u.= u.y * ( f - ( x * u.y * u.y ) );
228  u.= u.y * ( f - ( x * u.y * u.y ) );
229   
230  return number * u.y;
231}
232
233
234
235// Optimised squareroot computation (double version)
236double Xfast_squareroot_approximation(double number) {
237  double x;
238  const double f = 1.5;
239   
240  // Allow i and y to occupy the same memory   
241  union
242  {
243    long long i;
244    double y;
245  } u;
246
247  // Good initial guess
248  x = number * 0.5;
249  u.= number;
250  u.= 0x5fe6ec85e7de30daLL - ( u.i >> 1 );
251 
252  // Take a few iterations 
253  u.= u.y * ( f - ( x * u.y * u.y ) );
254  u.= u.y * ( f - ( x * u.y * u.y ) );
255
256  return number * u.y;
257}
258
259
260// Innermost flux function (using stage w=z+h)
261int _flux_function_central(double *q_left, double *q_right,
262                           double z_left, double z_right,
263                           double n1, double n2,
264                           double epsilon, 
265                           double h0,
266                           double limiting_threshold, 
267                           double g,
268                           double *edgeflux, double *max_speed) 
269{
270
271  /*Compute fluxes between volumes for the shallow water wave equation
272    cast in terms of the 'stage', w = h+z using
273    the 'central scheme' as described in
274
275    Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
276    Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
277    Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.
278
279    The implemented formula is given in equation (3.15) on page 714
280  */
281
282  int i;
283
284  double w_left, h_left, uh_left, vh_left, u_left;
285  double w_right, h_right, uh_right, vh_right, u_right;
286  double s_min, s_max, soundspeed_left, soundspeed_right;
287  double denom, inverse_denominator, z;
288
289  // Workspace (allocate once, use many)
290  static double q_left_rotated[3], q_right_rotated[3], flux_right[3], flux_left[3];
291
292  // Copy conserved quantities to protect from modification
293  q_left_rotated[0] = q_left[0];
294  q_right_rotated[0] = q_right[0];
295  q_left_rotated[1] = q_left[1];
296  q_right_rotated[1] = q_right[1];
297  q_left_rotated[2] = q_left[2];
298  q_right_rotated[2] = q_right[2];
299
300  // Align x- and y-momentum with x-axis
301  _rotate(q_left_rotated, n1, n2);
302  _rotate(q_right_rotated, n1, n2);
303
304  z = 0.5*(z_left + z_right); // Average elevation values.
305                              // Even though this will nominally allow
306                              // for discontinuities in the elevation data,
307                              // there is currently no numerical support for
308                              // this so results may be strange near
309                              // jumps in the bed.
310
311  // Compute speeds in x-direction
312  w_left = q_left_rotated[0];         
313  h_left = w_left - z;
314  uh_left = q_left_rotated[1];
315  u_left = _compute_speed(&uh_left, &h_left, 
316                          epsilon, h0, limiting_threshold);
317
318  w_right = q_right_rotated[0];
319  h_right = w_right - z;
320  uh_right = q_right_rotated[1];
321  u_right = _compute_speed(&uh_right, &h_right, 
322                           epsilon, h0, limiting_threshold);
323
324  // Momentum in y-direction
325  vh_left  = q_left_rotated[2];
326  vh_right = q_right_rotated[2];
327
328  // Limit y-momentum if necessary
329  // Leaving this out, improves speed significantly (Ole 27/5/2009)
330  // All validation tests pass, so do we really need it anymore?
331  _compute_speed(&vh_left, &h_left, 
332                 epsilon, h0, limiting_threshold);
333  _compute_speed(&vh_right, &h_right, 
334                 epsilon, h0, limiting_threshold);
335
336  // Maximal and minimal wave speeds
337  soundspeed_left  = sqrt(g*h_left);
338  soundspeed_right = sqrt(g*h_right); 
339 
340  // Code to use fast square root optimisation if desired.
341  // Timings on AMD 64 for the Okushiri profile gave the following timings
342  //
343  // SQRT           Total    Flux
344  //=============================
345  //
346  // Ref            405s     152s
347  // Fast (dbl)     453s     173s
348  // Fast (sng)     437s     171s
349  //
350  // Consequently, there is currently (14/5/2009) no reason to use this
351  // approximation.
352 
353  //soundspeed_left  = fast_squareroot_approximation(g*h_left);
354  //soundspeed_right = fast_squareroot_approximation(g*h_right);
355
356  s_max = max(u_left + soundspeed_left, u_right + soundspeed_right);
357  if (s_max < 0.0) 
358  {
359    s_max = 0.0;
360  }
361
362  s_min = min(u_left - soundspeed_left, u_right - soundspeed_right);
363  if (s_min > 0.0)
364  {
365    s_min = 0.0;
366  }
367 
368  // Flux formulas
369  flux_left[0] = u_left*h_left;
370  flux_left[1] = u_left*uh_left + 0.5*g*h_left*h_left;
371  flux_left[2] = u_left*vh_left;
372
373  flux_right[0] = u_right*h_right;
374  flux_right[1] = u_right*uh_right + 0.5*g*h_right*h_right;
375  flux_right[2] = u_right*vh_right;
376
377  // Flux computation
378  denom = s_max - s_min;
379  if (denom < epsilon) 
380  { // FIXME (Ole): Try using h0 here
381    memset(edgeflux, 0, 3*sizeof(double));
382    *max_speed = 0.0;
383  } 
384  else 
385  {
386    inverse_denominator = 1.0/denom;
387    for (i = 0; i < 3; i++) 
388    {
389      edgeflux[i] = s_max*flux_left[i] - s_min*flux_right[i];
390      edgeflux[i] += s_max*s_min*(q_right_rotated[i] - q_left_rotated[i]);
391      edgeflux[i] *= inverse_denominator;
392    }
393
394    // Maximal wavespeed
395    *max_speed = max(fabs(s_max), fabs(s_min));
396
397    // Rotate back
398    _rotate(edgeflux, n1, -n2);
399  }
400 
401  return 0;
402}
403
404double erfcc(double x){
405    double z,t,result;
406
407    z=fabs(x);
408    t=1.0/(1.0+0.5*z);
409    result=t*exp(-z*z-1.26551223+t*(1.00002368+t*(.37409196+
410         t*(.09678418+t*(-.18628806+t*(.27886807+t*(-1.13520398+
411         t*(1.48851587+t*(-.82215223+t*.17087277)))))))));
412    if (x < 0.0) result = 2.0-result;
413
414    return result;
415    }
416
417
418
419// Computational function for flux computation (using stage w=z+h)
420int flux_function_kinetic(double *q_left, double *q_right,
421          double z_left, double z_right,
422          double n1, double n2,
423          double epsilon, double H0, double g,
424          double *edgeflux, double *max_speed) {
425
426  /*Compute fluxes between volumes for the shallow water wave equation
427    cast in terms of the 'stage', w = h+z using
428    the 'central scheme' as described in
429
430    Zhang et. al., Advances in Water Resources, 26(6), 2003, 635-647.
431  */
432
433  int i;
434
435  double w_left, h_left, uh_left, vh_left, u_left, F_left;
436  double w_right, h_right, uh_right, vh_right, u_right, F_right;
437  double s_min, s_max, soundspeed_left, soundspeed_right;
438  double z;
439  double q_left_rotated[3], q_right_rotated[3];
440
441  double h0 = H0*H0; //This ensures a good balance when h approaches H0.
442  double limiting_threshold = 10*H0; // Avoid applying limiter below this 
443 
444  //Copy conserved quantities to protect from modification
445  for (i=0; i<3; i++) {
446    q_left_rotated[i] = q_left[i];
447    q_right_rotated[i] = q_right[i];
448  }
449
450  //Align x- and y-momentum with x-axis
451  _rotate(q_left_rotated, n1, n2);
452  _rotate(q_right_rotated, n1, n2);
453
454  z = (z_left+z_right)/2; //Take average of field values
455
456  //Compute speeds in x-direction
457  w_left = q_left_rotated[0];         
458  h_left = w_left-z;
459  uh_left = q_left_rotated[1];
460  u_left =_compute_speed(&uh_left, &h_left, 
461                         epsilon, h0, limiting_threshold);
462
463  w_right = q_right_rotated[0];
464  h_right = w_right-z;
465  uh_right = q_right_rotated[1];
466  u_right =_compute_speed(&uh_right, &h_right,
467                          epsilon, h0, limiting_threshold);
468
469
470  //Momentum in y-direction
471  vh_left  = q_left_rotated[2];
472  vh_right = q_right_rotated[2];
473
474
475  //Maximal and minimal wave speeds
476  soundspeed_left  = sqrt(g*h_left);
477  soundspeed_right = sqrt(g*h_right);
478
479  s_max = max(u_left+soundspeed_left, u_right+soundspeed_right);
480  if (s_max < 0.0) s_max = 0.0;
481
482  s_min = min(u_left-soundspeed_left, u_right-soundspeed_right);
483  if (s_min > 0.0) s_min = 0.0;
484
485
486  F_left  = 0.0;
487  F_right = 0.0;
488  if (h_left > 0.0) F_left = u_left/sqrt(g*h_left);
489  if (h_right > 0.0) F_right = u_right/sqrt(g*h_right);
490
491  for (i=0; i<3; i++) edgeflux[i] = 0.0;
492  *max_speed = 0.0;
493
494  edgeflux[0] = h_left*u_left/2.0*erfcc(-F_left) +  \
495          h_left*sqrt(g*h_left)/2.0/sqrt(pi)*exp(-(F_left*F_left)) + \
496          h_right*u_right/2.0*erfcc(F_right) -  \
497          h_right*sqrt(g*h_right)/2.0/sqrt(pi)*exp(-(F_right*F_right));
498
499  edgeflux[1] = (h_left*u_left*u_left + g/2.0*h_left*h_left)/2.0*erfcc(-F_left) + \
500          u_left*h_left*sqrt(g*h_left)/2.0/sqrt(pi)*exp(-(F_left*F_left)) + \
501          (h_right*u_right*u_right + g/2.0*h_right*h_right)/2.0*erfcc(F_right) -  \
502          u_right*h_right*sqrt(g*h_right)/2.0/sqrt(pi)*exp(-(F_right*F_right));
503
504  edgeflux[2] = vh_left*u_left/2.0*erfcc(-F_left) + \
505          vh_left*sqrt(g*h_left)/2.0/sqrt(pi)*exp(-(F_left*F_left)) + \
506          vh_right*u_right/2.0*erfcc(F_right) - \
507          vh_right*sqrt(g*h_right)/2.0/sqrt(pi)*exp(-(F_right*F_right));
508
509  //Maximal wavespeed
510  *max_speed = max(fabs(s_max), fabs(s_min));
511
512  //Rotate back
513  _rotate(edgeflux, n1, -n2);
514
515  return 0;
516}
517
518
519
520
521void _manning_friction_old(double g, double eps, int N,
522               double* w, double* zv,
523               double* uh, double* vh,
524               double* eta, double* xmom, double* ymom) {
525
526  int k, k3;
527  double S, h, z, z0, z1, z2;
528
529  for (k=0; k<N; k++) {
530    if (eta[k] > eps) {
531      k3 = 3*k;
532      // Get bathymetry
533      z0 = zv[k3 + 0];
534      z1 = zv[k3 + 1];
535      z2 = zv[k3 + 2];
536      z = (z0+z1+z2)/3.0;
537      h = w[k]-z;
538      if (h >= eps) {
539        S = -g * eta[k]*eta[k] * sqrt((uh[k]*uh[k] + vh[k]*vh[k]));
540        S /= pow(h, 7.0/3);      //Expensive (on Ole's home computer)
541        //S /= exp((7.0/3.0)*log(h));      //seems to save about 15% over manning_friction
542        //S /= h*h*(1 + h/3.0 - h*h/9.0); //FIXME: Could use a Taylor expansion
543
544
545        //Update momentum
546        xmom[k] += S*uh[k];
547        ymom[k] += S*vh[k];
548      }
549    }
550  }
551}
552
553
554void _manning_friction_new(double g, double eps, int N,
555                           double* x, double* w, double* zv,
556                           double* uh, double* vh,
557                           double* eta, double* xmom_update, double* ymom_update) {
558
559  int k, k3, k6;
560  double S, h, z, z0, z1, z2, zs, zx, zy;
561  double x0,y0,x1,y1,x2,y2;
562
563  for (k=0; k<N; k++) {
564    if (eta[k] > eps) {
565      k3 = 3*k;
566      // Get bathymetry
567      z0 = zv[k3 + 0];
568      z1 = zv[k3 + 1];
569      z2 = zv[k3 + 2];
570
571      // Compute bed slope
572      k6 = 6*k;  // base index
573 
574      x0 = x[k6 + 0];
575      y0 = x[k6 + 1];
576      x1 = x[k6 + 2];
577      y1 = x[k6 + 3];
578      x2 = x[k6 + 4];
579      y2 = x[k6 + 5];
580
581      _gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2, &zx, &zy);
582
583      zs = sqrt(1.0 + zx*zx + zy*zy);
584      z = (z0+z1+z2)/3.0;
585      h = w[k]-z;
586      if (h >= eps) {
587        S = -g * eta[k]*eta[k] * zs * sqrt((uh[k]*uh[k] + vh[k]*vh[k]));
588        S /= pow(h, 7.0/3);      //Expensive (on Ole's home computer)
589        //S /= exp((7.0/3.0)*log(h));      //seems to save about 15% over manning_friction
590        //S /= h*h*(1 + h/3.0 - h*h/9.0); //FIXME: Could use a Taylor expansion
591
592
593        //Update momentum
594        xmom_update[k] += S*uh[k];
595        ymom_update[k] += S*vh[k];
596      }
597    }
598  }
599}
600
601
602
603void _chezy_friction(double g, double eps, int N,
604                           double* x, double* w, double* zv,
605                           double* uh, double* vh,
606                           double* chezy, double* xmom_update, double* ymom_update) {
607
608  int k, k3, k6;
609  double S, h, z, z0, z1, z2, zs, zx, zy;
610  double x0,y0,x1,y1,x2,y2;
611
612  for (k=0; k<N; k++) {
613    if (chezy[k] > eps) {
614      k3 = 3*k;
615      // Get bathymetry
616      z0 = zv[k3 + 0];
617      z1 = zv[k3 + 1];
618      z2 = zv[k3 + 2];
619
620      // Compute bed slope
621      k6 = 6*k;  // base index
622
623      x0 = x[k6 + 0];
624      y0 = x[k6 + 1];
625      x1 = x[k6 + 2];
626      y1 = x[k6 + 3];
627      x2 = x[k6 + 4];
628      y2 = x[k6 + 5];
629
630      _gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2, &zx, &zy);
631
632      zs = sqrt(1.0 + zx*zx + zy*zy);
633      z = (z0+z1+z2)/3.0;
634      h = w[k]-z;
635      if (h >= eps) {
636        S = -g * chezy[k] * zs * sqrt((uh[k]*uh[k] + vh[k]*vh[k]));
637        S /= pow(h, 2.0);
638
639        //Update momentum
640        xmom_update[k] += S*uh[k];
641        ymom_update[k] += S*vh[k];
642      }
643    }
644  }
645}
646
647/*
648void _manning_friction_explicit(double g, double eps, int N,
649               double* w, double* z,
650               double* uh, double* vh,
651               double* eta, double* xmom, double* ymom) {
652
653  int k;
654  double S, h;
655
656  for (k=0; k<N; k++) {
657    if (eta[k] > eps) {
658      h = w[k]-z[k];
659      if (h >= eps) {
660    S = -g * eta[k]*eta[k] * sqrt((uh[k]*uh[k] + vh[k]*vh[k]));
661    S /= pow(h, 7.0/3);      //Expensive (on Ole's home computer)
662    //S /= exp(7.0/3.0*log(h));      //seems to save about 15% over manning_friction
663    //S /= h*h*(1 + h/3.0 - h*h/9.0); //FIXME: Could use a Taylor expansion
664
665
666    //Update momentum
667    xmom[k] += S*uh[k];
668    ymom[k] += S*vh[k];
669      }
670    }
671  }
672}
673*/
674
675
676
677
678double velocity_balance(double uh_i, 
679            double uh,
680            double h_i, 
681            double h, 
682            double alpha,
683            double epsilon) {
684  // Find alpha such that speed at the vertex is within one
685  // order of magnitude of the centroid speed   
686
687  // FIXME(Ole): Work in progress
688 
689  double a, b, estimate;
690  double m=10; // One order of magnitude - allow for velocity deviations at vertices
691
692 
693  printf("alpha = %f, uh_i=%f, uh=%f, h_i=%f, h=%f\n",
694     alpha, uh_i, uh, h_i, h);
695     
696 
697 
698   
699  // Shorthands and determine inequality
700  if (fabs(uh) < epsilon) {
701    a = 1.0e10; // Limit
702  } else {
703    a = fabs(uh_i - uh)/fabs(uh);
704  }
705     
706  if (h < epsilon) {
707    b = 1.0e10; // Limit
708  } else {   
709    b = m*fabs(h_i - h)/h;
710  }
711
712  printf("a %f, b %f\n", a, b); 
713
714  if (a > b) {
715    estimate = (m-1)/(a-b);     
716   
717    printf("Alpha %f, estimate %f\n", 
718       alpha, estimate);   
719       
720    if (alpha < estimate) {
721      printf("Adjusting alpha from %f to %f\n", 
722         alpha, estimate);
723      alpha = estimate;
724    }
725  } else {
726 
727    if (h < h_i) {
728      estimate = (m-1)/(a-b);             
729   
730      printf("Alpha %f, estimate %f\n", 
731         alpha, estimate);   
732       
733      if (alpha < estimate) {
734    printf("Adjusting alpha from %f to %f\n", 
735           alpha, estimate);
736    alpha = estimate;
737      }   
738    }
739    // Always fulfilled as alpha and m-1 are non negative
740  }
741 
742 
743  return alpha;
744}
745
746
747int _balance_deep_and_shallow(int N,
748                  double* wc,
749                  double* zc,
750                  double* wv,
751                  double* zv,
752                  double* hvbar, // Retire this
753                  double* xmomc,
754                  double* ymomc,
755                  double* xmomv,
756                  double* ymomv,
757                  double H0,
758                  int tight_slope_limiters,
759                  int use_centroid_velocities,
760                  double alpha_balance) {
761
762  int k, k3, i; 
763
764  double dz, hmin, alpha, h_diff, hc_k;
765  double epsilon = 1.0e-6; // FIXME: Temporary measure
766  double hv[3]; // Depths at vertices
767  double uc, vc; // Centroid speeds
768
769  // Compute linear combination between w-limited stages and
770  // h-limited stages close to the bed elevation.
771
772  for (k=0; k<N; k++) {
773    // Compute maximal variation in bed elevation
774    // This quantitiy is
775    //     dz = max_i abs(z_i - z_c)
776    // and it is independent of dimension
777    // In the 1d case zc = (z0+z1)/2
778    // In the 2d case zc = (z0+z1+z2)/3
779
780    k3 = 3*k;
781    hc_k = wc[k] - zc[k]; // Centroid value at triangle k
782   
783    dz = 0.0;
784    if (tight_slope_limiters == 0) {     
785      // FIXME: Try with this one precomputed
786      for (i=0; i<3; i++) {
787        dz = max(dz, fabs(zv[k3+i]-zc[k]));
788      }
789    }
790
791    // Calculate depth at vertices (possibly negative here!)
792    hv[0] = wv[k3] -   zv[k3];
793    hv[1] = wv[k3+1] - zv[k3+1];
794    hv[2] = wv[k3+2] - zv[k3+2];       
795   
796    // Calculate minimal depth across all three vertices
797    hmin = min(hv[0], min(hv[1], hv[2]));
798
799    //if (hmin < 0.0 ) {
800    //  printf("hmin = %f\n",hmin);
801    //}
802   
803
804    // Create alpha in [0,1], where alpha==0 means using the h-limited
805    // stage and alpha==1 means using the w-limited stage as
806    // computed by the gradient limiter (both 1st or 2nd order)
807    if (tight_slope_limiters == 0) {     
808      // If hmin > dz/alpha_balance then alpha = 1 and the bed will have no
809      // effect
810      // If hmin < 0 then alpha = 0 reverting to constant height above bed.
811      // The parameter alpha_balance==2 by default
812
813     
814      if (dz > 0.0) {
815        alpha = max( min( alpha_balance*hmin/dz, 1.0), 0.0 );     
816      } else {
817        alpha = 1.0;  // Flat bed
818      }
819      //printf("Using old style limiter\n");
820     
821    } else {
822
823      // Tight Slope Limiter (2007)
824   
825      // Make alpha as large as possible but still ensure that
826      // final depth is positive and that velocities at vertices
827      // are controlled
828   
829      if (hmin < H0) {
830        alpha = 1.0;
831        for (i=0; i<3; i++) {
832         
833          h_diff = hc_k - hv[i];     
834          if (h_diff <= 0) {
835            // Deep water triangle is further away from bed than
836            // shallow water (hbar < h). Any alpha will do
837     
838          } else { 
839            // Denominator is positive which means that we need some of the
840            // h-limited stage.
841           
842            alpha = min(alpha, (hc_k - H0)/h_diff);     
843          }
844        }
845
846        // Ensure alpha in [0,1]
847        if (alpha>1.0) alpha=1.0;
848        if (alpha<0.0) alpha=0.0;
849       
850      } else {
851        // Use w-limited stage exclusively in deeper water.
852        alpha = 1.0;       
853      }
854    }
855   
856   
857    //  Let
858    //
859    //    wvi be the w-limited stage (wvi = zvi + hvi)
860    //    wvi- be the h-limited state (wvi- = zvi + hvi-)
861    //
862    //
863    //  where i=0,1,2 denotes the vertex ids
864    //
865    //  Weighted balance between w-limited and h-limited stage is
866    //
867    //    wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi)
868    //
869    //  It follows that the updated wvi is
870    //    wvi := zvi + (1-alpha)*hvi- + alpha*hvi
871    //
872    //   Momentum is balanced between constant and limited
873
874   
875    if (alpha < 1) {     
876      for (i=0; i<3; i++) {
877       
878        wv[k3+i] = zv[k3+i] + (1-alpha)*hc_k + alpha*hv[i]; 
879
880        // Update momentum at vertices
881        if (use_centroid_velocities == 1) {
882          // This is a simple, efficient and robust option
883          // It uses first order approximation of velocities, but retains
884          // the order used by stage.
885   
886          // Speeds at centroids
887          if (hc_k > epsilon) {
888            uc = xmomc[k]/hc_k;
889            vc = ymomc[k]/hc_k;
890          } else {
891            uc = 0.0;
892            vc = 0.0;
893          }
894         
895          // Vertex momenta guaranteed to be consistent with depth guaranteeing
896          // controlled speed
897          hv[i] = wv[k3+i] - zv[k3+i]; // Recompute (balanced) vertex depth
898          xmomv[k3+i] = uc*hv[i];
899          ymomv[k3+i] = vc*hv[i];
900     
901        } else {
902          // Update momentum as a linear combination of
903          // xmomc and ymomc (shallow) and momentum
904          // from extrapolator xmomv and ymomv (deep).
905          // This assumes that values from xmomv and ymomv have
906          // been established e.g. by the gradient limiter.
907         
908          // FIXME (Ole): I think this should be used with vertex momenta
909          // computed above using centroid_velocities instead of xmomc
910          // and ymomc as they'll be more representative first order
911          // values.
912         
913          xmomv[k3+i] = (1-alpha)*xmomc[k] + alpha*xmomv[k3+i];
914          ymomv[k3+i] = (1-alpha)*ymomc[k] + alpha*ymomv[k3+i];
915         
916        }
917      }
918    }
919  }
920  return 0;
921}
922 
923
924
925
926int _protect(int N,
927         double minimum_allowed_height,
928         double maximum_allowed_speed,
929         double epsilon,
930         double* wc,
931         double* zc,
932         double* xmomc,
933         double* ymomc) {
934
935  int k;
936  double hc;
937  double u, v, reduced_speed;
938
939
940  // Protect against initesimal and negative heights 
941  if (maximum_allowed_speed < epsilon) {
942    for (k=0; k<N; k++) {
943      hc = wc[k] - zc[k];
944
945      if (hc < minimum_allowed_height) {
946       
947        // Set momentum to zero and ensure h is non negative
948        xmomc[k] = 0.0;
949        ymomc[k] = 0.0;
950        if (hc <= 0.0) wc[k] = zc[k];
951      }
952    }
953  } else {
954   
955    // Protect against initesimal and negative heights
956    for (k=0; k<N; k++) {
957      hc = wc[k] - zc[k];
958     
959      if (hc < minimum_allowed_height) {
960
961        //New code: Adjust momentum to guarantee speeds are physical
962        //          ensure h is non negative
963             
964        if (hc <= 0.0) {
965            wc[k] = zc[k];
966        xmomc[k] = 0.0;
967        ymomc[k] = 0.0;
968        } else {
969          //Reduce excessive speeds derived from division by small hc
970          //FIXME (Ole): This may be unnecessary with new slope limiters
971          //in effect.
972         
973          u = xmomc[k]/hc;
974          if (fabs(u) > maximum_allowed_speed) {
975            reduced_speed = maximum_allowed_speed * u/fabs(u);
976            //printf("Speed (u) has been reduced from %.3f to %.3f\n",
977            //   u, reduced_speed);
978            xmomc[k] = reduced_speed * hc;
979          }
980         
981          v = ymomc[k]/hc;
982          if (fabs(v) > maximum_allowed_speed) {
983            reduced_speed = maximum_allowed_speed * v/fabs(v);
984            //printf("Speed (v) has been reduced from %.3f to %.3f\n",
985            //   v, reduced_speed);
986            ymomc[k] = reduced_speed * hc;
987          }
988        }
989      }
990    }
991  }
992  return 0;
993}
994
995
996
997int _assign_wind_field_values(int N,
998                  double* xmom_update,
999                  double* ymom_update,
1000                  double* s_vec,
1001                  double* phi_vec,
1002                  double cw) {
1003
1004  // Assign windfield values to momentum updates
1005
1006  int k;
1007  double S, s, phi, u, v;
1008
1009  for (k=0; k<N; k++) {
1010
1011    s = s_vec[k];
1012    phi = phi_vec[k];
1013
1014    //Convert to radians
1015    phi = phi*pi/180;
1016
1017    //Compute velocity vector (u, v)
1018    u = s*cos(phi);
1019    v = s*sin(phi);
1020
1021    //Compute wind stress
1022    S = cw * sqrt(u*u + v*v);
1023    xmom_update[k] += S*u;
1024    ymom_update[k] += S*v;
1025  }
1026  return 0;
1027}
1028
1029
1030
1031///////////////////////////////////////////////////////////////////
1032// Gateways to Python
1033
1034
1035
1036PyObject *flux_function_central(PyObject *self, PyObject *args) {
1037  //
1038  // Gateway to innermost flux function.
1039  // This is only used by the unit tests as the c implementation is
1040  // normally called by compute_fluxes in this module.
1041
1042
1043  PyArrayObject *normal, *ql, *qr,  *edgeflux;
1044  double g, epsilon, max_speed, H0, zl, zr;
1045  double h0, limiting_threshold;
1046
1047  if (!PyArg_ParseTuple(args, "OOOddOddd",
1048            &normal, &ql, &qr, &zl, &zr, &edgeflux,
1049            &epsilon, &g, &H0)) {
1050      report_python_error(AT, "could not parse input arguments");
1051      return NULL;
1052  }
1053
1054 
1055  h0 = H0*H0; // This ensures a good balance when h approaches H0.
1056              // But evidence suggests that h0 can be as little as
1057              // epsilon!
1058             
1059  limiting_threshold = 10*H0; // Avoid applying limiter below this
1060                              // threshold for performance reasons.
1061                              // See ANUGA manual under flux limiting 
1062 
1063  _flux_function_central((double*) ql -> data, 
1064                         (double*) qr -> data, 
1065                         zl, 
1066                         zr,                         
1067                         ((double*) normal -> data)[0],
1068                         ((double*) normal -> data)[1],         
1069                         epsilon, h0, limiting_threshold,
1070                         g,
1071                         (double*) edgeflux -> data, 
1072                         &max_speed);
1073 
1074  return Py_BuildValue("d", max_speed); 
1075}
1076
1077
1078
1079PyObject *gravity(PyObject *self, PyObject *args) {
1080  //
1081  //  gravity(g, h, v, x, xmom, ymom)
1082  //
1083
1084
1085  PyArrayObject *h, *z, *x, *xmom, *ymom;
1086  int k, N, k3, k6;
1087  double g, avg_h, zx, zy;
1088  double x0, y0, x1, y1, x2, y2, z0, z1, z2;
1089  //double epsilon;
1090
1091  if (!PyArg_ParseTuple(args, "dOOOOO",
1092            &g, &h, &z, &x,
1093            &xmom, &ymom)) {
1094    //&epsilon)) {
1095      report_python_error(AT, "could not parse input arguments");
1096      return NULL;
1097  }
1098
1099  // check that numpy array objects arrays are C contiguous memory
1100  CHECK_C_CONTIG(h);
1101  CHECK_C_CONTIG(z);
1102  CHECK_C_CONTIG(x);
1103  CHECK_C_CONTIG(xmom);
1104  CHECK_C_CONTIG(ymom);
1105
1106  N = h -> dimensions[0];
1107  for (k=0; k<N; k++) {
1108    k3 = 3*k;  // base index
1109
1110    // Get bathymetry
1111    z0 = ((double*) z -> data)[k3 + 0];
1112    z1 = ((double*) z -> data)[k3 + 1];
1113    z2 = ((double*) z -> data)[k3 + 2];
1114
1115    // Optimise for flat bed
1116    // Note (Ole): This didn't produce measurable speed up.
1117    // Revisit later
1118    //if (fabs(z0-z1)<epsilon && fabs(z1-z2)<epsilon) {
1119    //  continue;
1120    //}
1121
1122    // Get average depth from centroid values
1123    avg_h = ((double *) h -> data)[k];
1124
1125    // Compute bed slope
1126    k6 = 6*k;  // base index
1127 
1128    x0 = ((double*) x -> data)[k6 + 0];
1129    y0 = ((double*) x -> data)[k6 + 1];
1130    x1 = ((double*) x -> data)[k6 + 2];
1131    y1 = ((double*) x -> data)[k6 + 3];
1132    x2 = ((double*) x -> data)[k6 + 4];
1133    y2 = ((double*) x -> data)[k6 + 5];
1134
1135
1136    _gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2, &zx, &zy);
1137
1138    // Update momentum
1139    ((double*) xmom -> data)[k] += -g*zx*avg_h;
1140    ((double*) ymom -> data)[k] += -g*zy*avg_h;
1141  }
1142
1143  return Py_BuildValue("");
1144}
1145
1146
1147PyObject *manning_friction_new(PyObject *self, PyObject *args) {
1148  //
1149  // manning_friction_new(g, eps, x, h, uh, vh, z, eta, xmom_update, ymom_update)
1150  //
1151
1152
1153  PyArrayObject *x, *w, *z, *uh, *vh, *eta, *xmom, *ymom;
1154  int N;
1155  double g, eps;
1156
1157  if (!PyArg_ParseTuple(args, "ddOOOOOOOO",
1158                        &g, &eps, &x, &w, &uh, &vh, &z,  &eta, &xmom, &ymom)) {
1159      report_python_error(AT, "could not parse input arguments");
1160      return NULL;
1161  }
1162
1163  // check that numpy array objects arrays are C contiguous memory
1164  CHECK_C_CONTIG(x);
1165  CHECK_C_CONTIG(w);
1166  CHECK_C_CONTIG(z);
1167  CHECK_C_CONTIG(uh);
1168  CHECK_C_CONTIG(vh);
1169  CHECK_C_CONTIG(eta);
1170  CHECK_C_CONTIG(xmom);
1171  CHECK_C_CONTIG(ymom);
1172
1173  N = w -> dimensions[0];
1174
1175  _manning_friction_new(g, eps, N,
1176                        (double*) x -> data,
1177                        (double*) w -> data,
1178                        (double*) z -> data,
1179                        (double*) uh -> data,
1180                        (double*) vh -> data,
1181                        (double*) eta -> data,
1182                        (double*) xmom -> data,
1183                        (double*) ymom -> data);
1184
1185  return Py_BuildValue("");
1186}
1187
1188
1189PyObject *chezy_friction(PyObject *self, PyObject *args) {
1190  //
1191  // chezy_friction(g, eps, x, h, uh, vh, z, chezy, xmom_update, ymom_update)
1192  //
1193
1194
1195  PyArrayObject *x, *w, *z, *uh, *vh, *chezy, *xmom, *ymom;
1196  int N;
1197  double g, eps;
1198
1199  if (!PyArg_ParseTuple(args, "ddOOOOOOOO",
1200                        &g, &eps, &x, &w, &uh, &vh, &z,  &chezy, &xmom, &ymom)) {
1201      report_python_error(AT, "could not parse input arguments");
1202      return NULL;
1203  }
1204
1205  // check that numpy array objects arrays are C contiguous memory
1206  CHECK_C_CONTIG(x);
1207  CHECK_C_CONTIG(w);
1208  CHECK_C_CONTIG(z);
1209  CHECK_C_CONTIG(uh);
1210  CHECK_C_CONTIG(vh);
1211  CHECK_C_CONTIG(chezy);
1212  CHECK_C_CONTIG(xmom);
1213  CHECK_C_CONTIG(ymom);
1214
1215  N = w -> dimensions[0];
1216
1217  _chezy_friction(g, eps, N,
1218                        (double*) x -> data,
1219                        (double*) w -> data,
1220                        (double*) z -> data,
1221                        (double*) uh -> data,
1222                        (double*) vh -> data,
1223                        (double*) chezy -> data,
1224                        (double*) xmom -> data,
1225                        (double*) ymom -> data);
1226
1227  return Py_BuildValue("");
1228}
1229
1230
1231PyObject *manning_friction_old(PyObject *self, PyObject *args) {
1232  //
1233  // manning_friction_old(g, eps, h, uh, vh, z, eta, xmom_update, ymom_update)
1234  //
1235
1236
1237  PyArrayObject *w, *z, *uh, *vh, *eta, *xmom, *ymom;
1238  int N;
1239  double g, eps;
1240
1241  if (!PyArg_ParseTuple(args, "ddOOOOOOO",
1242                        &g, &eps, &w, &uh, &vh, &z,  &eta, &xmom, &ymom)) {
1243      report_python_error(AT, "could not parse input arguments");
1244      return NULL;
1245  }
1246
1247  // check that numpy array objects arrays are C contiguous memory
1248  CHECK_C_CONTIG(w);
1249  CHECK_C_CONTIG(z);
1250  CHECK_C_CONTIG(uh);
1251  CHECK_C_CONTIG(vh);
1252  CHECK_C_CONTIG(eta);
1253  CHECK_C_CONTIG(xmom);
1254  CHECK_C_CONTIG(ymom);
1255
1256  N = w -> dimensions[0];
1257
1258  _manning_friction_old(g, eps, N,
1259                        (double*) w -> data,
1260                        (double*) z -> data,
1261                        (double*) uh -> data,
1262                        (double*) vh -> data,
1263                        (double*) eta -> data,
1264                        (double*) xmom -> data,
1265                        (double*) ymom -> data);
1266
1267  return Py_BuildValue("");
1268}
1269
1270
1271/*
1272PyObject *manning_friction_explicit(PyObject *self, PyObject *args) {
1273  //
1274  // manning_friction_explicit(g, eps, h, uh, vh, eta, xmom_update, ymom_update)
1275  //
1276
1277
1278  PyArrayObject *w, *z, *uh, *vh, *eta, *xmom, *ymom;
1279  int N;
1280  double g, eps;
1281
1282  if (!PyArg_ParseTuple(args, "ddOOOOOOO",
1283            &g, &eps, &w, &z, &uh, &vh, &eta,
1284            &xmom, &ymom))
1285    return NULL;
1286
1287  // check that numpy array objects arrays are C contiguous memory
1288  CHECK_C_CONTIG(w);
1289  CHECK_C_CONTIG(z);
1290  CHECK_C_CONTIG(uh);
1291  CHECK_C_CONTIG(vh);
1292  CHECK_C_CONTIG(eta);
1293  CHECK_C_CONTIG(xmom);
1294  CHECK_C_CONTIG(ymom);
1295
1296  N = w -> dimensions[0];
1297  _manning_friction_explicit(g, eps, N,
1298            (double*) w -> data,
1299            (double*) z -> data,
1300            (double*) uh -> data,
1301            (double*) vh -> data,
1302            (double*) eta -> data,
1303            (double*) xmom -> data,
1304            (double*) ymom -> data);
1305
1306  return Py_BuildValue("");
1307}
1308*/
1309
1310
1311
1312// Computational routine
1313int _extrapolate_second_order_sw(int number_of_elements,
1314                                 double epsilon,
1315                                 double minimum_allowed_height,
1316                                 double beta_w,
1317                                 double beta_w_dry,
1318                                 double beta_uh,
1319                                 double beta_uh_dry,
1320                                 double beta_vh,
1321                                 double beta_vh_dry,
1322                                 long* surrogate_neighbours,
1323                                 long* number_of_boundaries,
1324                                 double* centroid_coordinates,
1325                                 double* stage_centroid_values,
1326                                 double* xmom_centroid_values,
1327                                 double* ymom_centroid_values,
1328                                 double* elevation_centroid_values,
1329                                 double* vertex_coordinates,
1330                                 double* stage_vertex_values,
1331                                 double* xmom_vertex_values,
1332                                 double* ymom_vertex_values,
1333                                 double* elevation_vertex_values,
1334                                 int optimise_dry_cells) {
1335                 
1336                 
1337
1338  // Local variables
1339  double a, b; // Gradient vector used to calculate vertex values from centroids
1340  int k, k0, k1, k2, k3, k6, coord_index, i;
1341  double x, y, x0, y0, x1, y1, x2, y2, xv0, yv0, xv1, yv1, xv2, yv2; // Vertices of the auxiliary triangle
1342  double dx1, dx2, dy1, dy2, dxv0, dxv1, dxv2, dyv0, dyv1, dyv2, dq0, dq1, dq2, area2, inv_area2;
1343  double dqv[3], qmin, qmax, hmin, hmax;
1344  double hc, h0, h1, h2, beta_tmp, hfactor;
1345
1346   
1347  for (k = 0; k < number_of_elements; k++) 
1348  {
1349    k3=k*3;
1350    k6=k*6;
1351
1352    if (number_of_boundaries[k]==3)
1353    {
1354      // No neighbours, set gradient on the triangle to zero
1355     
1356      stage_vertex_values[k3]   = stage_centroid_values[k];
1357      stage_vertex_values[k3+1] = stage_centroid_values[k];
1358      stage_vertex_values[k3+2] = stage_centroid_values[k];
1359      xmom_vertex_values[k3]    = xmom_centroid_values[k];
1360      xmom_vertex_values[k3+1]  = xmom_centroid_values[k];
1361      xmom_vertex_values[k3+2]  = xmom_centroid_values[k];
1362      ymom_vertex_values[k3]    = ymom_centroid_values[k];
1363      ymom_vertex_values[k3+1]  = ymom_centroid_values[k];
1364      ymom_vertex_values[k3+2]  = ymom_centroid_values[k];
1365     
1366      continue;
1367    }
1368    else 
1369    {
1370      // Triangle k has one or more neighbours.
1371      // Get centroid and vertex coordinates of the triangle
1372     
1373      // Get the vertex coordinates
1374      xv0 = vertex_coordinates[k6];   
1375      yv0 = vertex_coordinates[k6+1];
1376      xv1 = vertex_coordinates[k6+2]; 
1377      yv1 = vertex_coordinates[k6+3];
1378      xv2 = vertex_coordinates[k6+4]; 
1379      yv2 = vertex_coordinates[k6+5];
1380     
1381      // Get the centroid coordinates
1382      coord_index = 2*k;
1383      x = centroid_coordinates[coord_index];
1384      y = centroid_coordinates[coord_index+1];
1385     
1386      // Store x- and y- differentials for the vertices of
1387      // triangle k relative to the centroid
1388      dxv0 = xv0 - x; 
1389      dxv1 = xv1 - x; 
1390      dxv2 = xv2 - x;
1391      dyv0 = yv0 - y; 
1392      dyv1 = yv1 - y; 
1393      dyv2 = yv2 - y;
1394    }
1395
1396
1397           
1398    if (number_of_boundaries[k]<=1)
1399    {
1400      //==============================================
1401      // Number of boundaries <= 1
1402      //==============================================   
1403   
1404   
1405      // If no boundaries, auxiliary triangle is formed
1406      // from the centroids of the three neighbours
1407      // If one boundary, auxiliary triangle is formed
1408      // from this centroid and its two neighbours
1409     
1410      k0 = surrogate_neighbours[k3];
1411      k1 = surrogate_neighbours[k3 + 1];
1412      k2 = surrogate_neighbours[k3 + 2];
1413     
1414      // Get the auxiliary triangle's vertex coordinates
1415      // (really the centroids of neighbouring triangles)
1416      coord_index = 2*k0;
1417      x0 = centroid_coordinates[coord_index];
1418      y0 = centroid_coordinates[coord_index+1];
1419     
1420      coord_index = 2*k1;
1421      x1 = centroid_coordinates[coord_index];
1422      y1 = centroid_coordinates[coord_index+1];
1423     
1424      coord_index = 2*k2;
1425      x2 = centroid_coordinates[coord_index];
1426      y2 = centroid_coordinates[coord_index+1];
1427     
1428      // Store x- and y- differentials for the vertices
1429      // of the auxiliary triangle
1430      dx1 = x1 - x0; 
1431      dx2 = x2 - x0;
1432      dy1 = y1 - y0; 
1433      dy2 = y2 - y0;
1434     
1435      // Calculate 2*area of the auxiliary triangle
1436      // The triangle is guaranteed to be counter-clockwise     
1437      area2 = dy2*dx1 - dy1*dx2; 
1438     
1439      // If the mesh is 'weird' near the boundary,
1440      // the triangle might be flat or clockwise:
1441      if (area2 <= 0) 
1442      {
1443          report_python_error(AT, "Negative triangle area");
1444          return -1;
1445      } 
1446     
1447      // Calculate heights of neighbouring cells
1448      hc = stage_centroid_values[k]  - elevation_centroid_values[k];
1449      h0 = stage_centroid_values[k0] - elevation_centroid_values[k0];
1450      h1 = stage_centroid_values[k1] - elevation_centroid_values[k1];
1451      h2 = stage_centroid_values[k2] - elevation_centroid_values[k2];
1452      hmin = min(min(h0, min(h1, h2)), hc);
1453      //hfactor = hc/(hc + 1.0);
1454
1455      hfactor = 0.0;
1456      if (hmin > 0.001) 
1457      {
1458        hfactor = (hmin - 0.001)/(hmin + 0.004);
1459      }
1460     
1461      if (optimise_dry_cells) 
1462      {     
1463        // Check if linear reconstruction is necessary for triangle k
1464        // This check will exclude dry cells.
1465
1466        hmax = max(h0, max(h1, h2));     
1467        if (hmax < epsilon) 
1468        {
1469            continue;
1470        }
1471      }
1472   
1473      //-----------------------------------
1474      // stage
1475      //-----------------------------------     
1476     
1477      // Calculate the difference between vertex 0 of the auxiliary
1478      // triangle and the centroid of triangle k
1479      dq0 = stage_centroid_values[k0] - stage_centroid_values[k];
1480     
1481      // Calculate differentials between the vertices
1482      // of the auxiliary triangle (centroids of neighbouring triangles)
1483      dq1 = stage_centroid_values[k1] - stage_centroid_values[k0];
1484      dq2 = stage_centroid_values[k2] - stage_centroid_values[k0];
1485     
1486          inv_area2 = 1.0/area2;
1487      // Calculate the gradient of stage on the auxiliary triangle
1488      a = dy2*dq1 - dy1*dq2;
1489      a *= inv_area2;
1490      b = dx1*dq2 - dx2*dq1;
1491      b *= inv_area2;
1492     
1493      // Calculate provisional jumps in stage from the centroid
1494      // of triangle k to its vertices, to be limited
1495      dqv[0] = a*dxv0 + b*dyv0;
1496      dqv[1] = a*dxv1 + b*dyv1;
1497      dqv[2] = a*dxv2 + b*dyv2;
1498     
1499      // Now we want to find min and max of the centroid and the
1500      // vertices of the auxiliary triangle and compute jumps
1501      // from the centroid to the min and max
1502      find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
1503     
1504      // Playing with dry wet interface
1505      //hmin = qmin;
1506      //beta_tmp = beta_w_dry;
1507      //if (hmin>minimum_allowed_height)
1508      beta_tmp = beta_w_dry + (beta_w - beta_w_dry) * hfactor;
1509   
1510      //printf("min_alled_height = %f\n",minimum_allowed_height);
1511      //printf("hmin = %f\n",hmin);
1512      //printf("beta_w = %f\n",beta_w);
1513      //printf("beta_tmp = %f\n",beta_tmp);
1514      // Limit the gradient
1515      limit_gradient(dqv, qmin, qmax, beta_tmp); 
1516     
1517      //for (i=0;i<3;i++)
1518      stage_vertex_values[k3+0] = stage_centroid_values[k] + dqv[0];
1519          stage_vertex_values[k3+1] = stage_centroid_values[k] + dqv[1];
1520          stage_vertex_values[k3+2] = stage_centroid_values[k] + dqv[2];
1521     
1522     
1523      //-----------------------------------
1524      // xmomentum
1525      //-----------------------------------           
1526
1527      // Calculate the difference between vertex 0 of the auxiliary
1528      // triangle and the centroid of triangle k     
1529      dq0 = xmom_centroid_values[k0] - xmom_centroid_values[k];
1530     
1531      // Calculate differentials between the vertices
1532      // of the auxiliary triangle
1533      dq1 = xmom_centroid_values[k1] - xmom_centroid_values[k0];
1534      dq2 = xmom_centroid_values[k2] - xmom_centroid_values[k0];
1535     
1536      // Calculate the gradient of xmom on the auxiliary triangle
1537      a = dy2*dq1 - dy1*dq2;
1538      a *= inv_area2;
1539      b = dx1*dq2 - dx2*dq1;
1540      b *= inv_area2;
1541     
1542      // Calculate provisional jumps in stage from the centroid
1543      // of triangle k to its vertices, to be limited     
1544      dqv[0] = a*dxv0+b*dyv0;
1545      dqv[1] = a*dxv1+b*dyv1;
1546      dqv[2] = a*dxv2+b*dyv2;
1547     
1548      // Now we want to find min and max of the centroid and the
1549      // vertices of the auxiliary triangle and compute jumps
1550      // from the centroid to the min and max
1551      find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
1552      //beta_tmp = beta_uh;
1553      //if (hmin<minimum_allowed_height)
1554      //beta_tmp = beta_uh_dry;
1555      beta_tmp = beta_uh_dry + (beta_uh - beta_uh_dry) * hfactor;
1556
1557      // Limit the gradient
1558      limit_gradient(dqv, qmin, qmax, beta_tmp);
1559
1560      for (i=0; i < 3; i++)
1561      {
1562        xmom_vertex_values[k3+i] = xmom_centroid_values[k] + dqv[i];
1563      }
1564     
1565      //-----------------------------------
1566      // ymomentum
1567      //-----------------------------------                 
1568
1569      // Calculate the difference between vertex 0 of the auxiliary
1570      // triangle and the centroid of triangle k
1571      dq0 = ymom_centroid_values[k0] - ymom_centroid_values[k];
1572     
1573      // Calculate differentials between the vertices
1574      // of the auxiliary triangle
1575      dq1 = ymom_centroid_values[k1] - ymom_centroid_values[k0];
1576      dq2 = ymom_centroid_values[k2] - ymom_centroid_values[k0];
1577     
1578      // Calculate the gradient of xmom on the auxiliary triangle
1579      a = dy2*dq1 - dy1*dq2;
1580      a *= inv_area2;
1581      b = dx1*dq2 - dx2*dq1;
1582      b *= inv_area2;
1583     
1584      // Calculate provisional jumps in stage from the centroid
1585      // of triangle k to its vertices, to be limited
1586      dqv[0] = a*dxv0 + b*dyv0;
1587      dqv[1] = a*dxv1 + b*dyv1;
1588      dqv[2] = a*dxv2 + b*dyv2;
1589     
1590      // Now we want to find min and max of the centroid and the
1591      // vertices of the auxiliary triangle and compute jumps
1592      // from the centroid to the min and max
1593      find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
1594     
1595      //beta_tmp = beta_vh;
1596      //
1597      //if (hmin<minimum_allowed_height)
1598      //beta_tmp = beta_vh_dry;
1599      beta_tmp = beta_vh_dry + (beta_vh - beta_vh_dry) * hfactor;   
1600
1601      // Limit the gradient
1602      limit_gradient(dqv, qmin, qmax, beta_tmp);
1603     
1604      for (i=0;i<3;i++)
1605      {
1606        ymom_vertex_values[k3 + i] = ymom_centroid_values[k] + dqv[i];
1607      }
1608    } // End number_of_boundaries <=1
1609    else
1610    {
1611
1612      //==============================================
1613      // Number of boundaries == 2
1614      //==============================================       
1615       
1616      // One internal neighbour and gradient is in direction of the neighbour's centroid
1617     
1618      // Find the only internal neighbour (k1?)
1619      for (k2 = k3; k2 < k3 + 3; k2++)
1620      {
1621      // Find internal neighbour of triangle k     
1622      // k2 indexes the edges of triangle k   
1623     
1624          if (surrogate_neighbours[k2] != k)
1625          {
1626             break;
1627          }
1628      }
1629     
1630      if ((k2 == k3 + 3)) 
1631      {
1632        // If we didn't find an internal neighbour
1633        report_python_error(AT, "Internal neighbour not found");
1634        return -1;
1635      }
1636     
1637      k1 = surrogate_neighbours[k2];
1638     
1639      // The coordinates of the triangle are already (x,y).
1640      // Get centroid of the neighbour (x1,y1)
1641      coord_index = 2*k1;
1642      x1 = centroid_coordinates[coord_index];
1643      y1 = centroid_coordinates[coord_index + 1];
1644     
1645      // Compute x- and y- distances between the centroid of
1646      // triangle k and that of its neighbour
1647      dx1 = x1 - x; 
1648      dy1 = y1 - y;
1649     
1650      // Set area2 as the square of the distance
1651      area2 = dx1*dx1 + dy1*dy1;
1652     
1653      // Set dx2=(x1-x0)/((x1-x0)^2+(y1-y0)^2)
1654      // and dy2=(y1-y0)/((x1-x0)^2+(y1-y0)^2) which
1655      // respectively correspond to the x- and y- gradients
1656      // of the conserved quantities
1657      dx2 = 1.0/area2;
1658      dy2 = dx2*dy1;
1659      dx2 *= dx1;
1660     
1661     
1662      //-----------------------------------
1663      // stage
1664      //-----------------------------------           
1665
1666      // Compute differentials
1667      dq1 = stage_centroid_values[k1] - stage_centroid_values[k];
1668     
1669      // Calculate the gradient between the centroid of triangle k
1670      // and that of its neighbour
1671      a = dq1*dx2;
1672      b = dq1*dy2;
1673     
1674      // Calculate provisional vertex jumps, to be limited
1675      dqv[0] = a*dxv0 + b*dyv0;
1676      dqv[1] = a*dxv1 + b*dyv1;
1677      dqv[2] = a*dxv2 + b*dyv2;
1678     
1679      // Now limit the jumps
1680      if (dq1>=0.0)
1681      {
1682        qmin=0.0;
1683        qmax=dq1;
1684      }
1685      else
1686      {
1687        qmin = dq1;
1688        qmax = 0.0;
1689      }
1690     
1691      // Limit the gradient
1692      limit_gradient(dqv, qmin, qmax, beta_w);
1693     
1694      //for (i=0; i < 3; i++)
1695      //{
1696      stage_vertex_values[k3] = stage_centroid_values[k] + dqv[0];
1697      stage_vertex_values[k3 + 1] = stage_centroid_values[k] + dqv[1];
1698      stage_vertex_values[k3 + 2] = stage_centroid_values[k] + dqv[2];
1699      //}
1700     
1701      //-----------------------------------
1702      // xmomentum
1703      //-----------------------------------                       
1704     
1705      // Compute differentials
1706      dq1 = xmom_centroid_values[k1] - xmom_centroid_values[k];
1707     
1708      // Calculate the gradient between the centroid of triangle k
1709      // and that of its neighbour
1710      a = dq1*dx2;
1711      b = dq1*dy2;
1712     
1713      // Calculate provisional vertex jumps, to be limited
1714      dqv[0] = a*dxv0+b*dyv0;
1715      dqv[1] = a*dxv1+b*dyv1;
1716      dqv[2] = a*dxv2+b*dyv2;
1717     
1718      // Now limit the jumps
1719      if (dq1 >= 0.0)
1720      {
1721        qmin = 0.0;
1722        qmax = dq1;
1723      }
1724      else
1725      {
1726        qmin = dq1;
1727        qmax = 0.0;
1728      }
1729     
1730      // Limit the gradient     
1731      limit_gradient(dqv, qmin, qmax, beta_w);
1732     
1733      //for (i=0;i<3;i++)
1734      //xmom_vertex_values[k3] = xmom_centroid_values[k] + dqv[0];
1735      //xmom_vertex_values[k3 + 1] = xmom_centroid_values[k] + dqv[1];
1736      //xmom_vertex_values[k3 + 2] = xmom_centroid_values[k] + dqv[2];
1737     
1738      for (i = 0; i < 3;i++)
1739      {
1740          xmom_vertex_values[k3 + i] = xmom_centroid_values[k] + dqv[i];
1741      }
1742     
1743      //-----------------------------------
1744      // ymomentum
1745      //-----------------------------------                       
1746
1747      // Compute differentials
1748      dq1 = ymom_centroid_values[k1] - ymom_centroid_values[k];
1749     
1750      // Calculate the gradient between the centroid of triangle k
1751      // and that of its neighbour
1752      a = dq1*dx2;
1753      b = dq1*dy2;
1754     
1755      // Calculate provisional vertex jumps, to be limited
1756      dqv[0] = a*dxv0 + b*dyv0;
1757      dqv[1] = a*dxv1 + b*dyv1;
1758      dqv[2] = a*dxv2 + b*dyv2;
1759     
1760      // Now limit the jumps
1761      if (dq1>=0.0)
1762      {
1763        qmin = 0.0;
1764        qmax = dq1;
1765      } 
1766      else
1767      {
1768        qmin = dq1;
1769        qmax = 0.0;
1770      }
1771     
1772      // Limit the gradient
1773      limit_gradient(dqv, qmin, qmax, beta_w);
1774     
1775      //for (i=0;i<3;i++)
1776      //ymom_vertex_values[k3] = ymom_centroid_values[k] + dqv[0];
1777      //ymom_vertex_values[k3 + 1] = ymom_centroid_values[k] + dqv[1];
1778      //ymom_vertex_values[k3 + 2] = ymom_centroid_values[k] + dqv[2];
1779
1780for (i=0;i<3;i++)
1781        {
1782ymom_vertex_values[k3 + i] = ymom_centroid_values[k] + dqv[i];
1783        }
1784//ymom_vertex_values[k3] = ymom_centroid_values[k] + dqv[0];
1785//ymom_vertex_values[k3 + 1] = ymom_centroid_values[k] + dqv[1];
1786//ymom_vertex_values[k3 + 2] = ymom_centroid_values[k] + dqv[2];
1787    } // else [number_of_boundaries==2]
1788  } // for k=0 to number_of_elements-1
1789 
1790  return 0;
1791}           
1792
1793
1794PyObject *extrapolate_second_order_sw(PyObject *self, PyObject *args) {
1795  /*Compute the vertex values based on a linear reconstruction
1796    on each triangle
1797   
1798    These values are calculated as follows:
1799    1) For each triangle not adjacent to a boundary, we consider the
1800       auxiliary triangle formed by the centroids of its three
1801       neighbours.
1802    2) For each conserved quantity, we integrate around the auxiliary
1803       triangle's boundary the product of the quantity and the outward
1804       normal vector. Dividing by the triangle area gives (a,b), the
1805       average of the vector (q_x,q_y) on the auxiliary triangle.
1806       We suppose that the linear reconstruction on the original
1807       triangle has gradient (a,b).
1808    3) Provisional vertex jumps dqv[0,1,2] are computed and these are
1809       then limited by calling the functions find_qmin_and_qmax and
1810       limit_gradient
1811
1812    Python call:
1813    extrapolate_second_order_sw(domain.surrogate_neighbours,
1814                                domain.number_of_boundaries
1815                                domain.centroid_coordinates,
1816                                Stage.centroid_values
1817                                Xmom.centroid_values
1818                                Ymom.centroid_values
1819                                domain.vertex_coordinates,
1820                                Stage.vertex_values,
1821                                Xmom.vertex_values,
1822                                Ymom.vertex_values)
1823
1824    Post conditions:
1825            The vertices of each triangle have values from a
1826        limited linear reconstruction
1827        based on centroid values
1828
1829  */
1830  PyArrayObject *surrogate_neighbours,
1831    *number_of_boundaries,
1832    *centroid_coordinates,
1833    *stage_centroid_values,
1834    *xmom_centroid_values,
1835    *ymom_centroid_values,
1836    *elevation_centroid_values,
1837    *vertex_coordinates,
1838    *stage_vertex_values,
1839    *xmom_vertex_values,
1840    *ymom_vertex_values,
1841    *elevation_vertex_values;
1842 
1843  PyObject *domain;
1844
1845 
1846  double beta_w, beta_w_dry, beta_uh, beta_uh_dry, beta_vh, beta_vh_dry;   
1847  double minimum_allowed_height, epsilon;
1848  int optimise_dry_cells, number_of_elements, e;
1849 
1850  // Provisional jumps from centroids to v'tices and safety factor re limiting
1851  // by which these jumps are limited
1852  // Convert Python arguments to C
1853  if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOi",
1854            &domain,
1855            &surrogate_neighbours,
1856            &number_of_boundaries,
1857            &centroid_coordinates,
1858            &stage_centroid_values,
1859            &xmom_centroid_values,
1860            &ymom_centroid_values,
1861            &elevation_centroid_values,
1862            &vertex_coordinates,
1863            &stage_vertex_values,
1864            &xmom_vertex_values,
1865            &ymom_vertex_values,
1866            &elevation_vertex_values,
1867            &optimise_dry_cells)) {         
1868
1869    report_python_error(AT, "could not parse input arguments");
1870    return NULL;
1871  }
1872
1873  // check that numpy array objects arrays are C contiguous memory
1874  CHECK_C_CONTIG(surrogate_neighbours);
1875  CHECK_C_CONTIG(number_of_boundaries);
1876  CHECK_C_CONTIG(centroid_coordinates);
1877  CHECK_C_CONTIG(stage_centroid_values);
1878  CHECK_C_CONTIG(xmom_centroid_values);
1879  CHECK_C_CONTIG(ymom_centroid_values);
1880  CHECK_C_CONTIG(elevation_centroid_values);
1881  CHECK_C_CONTIG(vertex_coordinates);
1882  CHECK_C_CONTIG(stage_vertex_values);
1883  CHECK_C_CONTIG(xmom_vertex_values);
1884  CHECK_C_CONTIG(ymom_vertex_values);
1885  CHECK_C_CONTIG(elevation_vertex_values);
1886 
1887  // Get the safety factor beta_w, set in the config.py file.
1888  // This is used in the limiting process
1889 
1890
1891  beta_w                 = get_python_double(domain,"beta_w");
1892  beta_w_dry             = get_python_double(domain,"beta_w_dry");
1893  beta_uh                = get_python_double(domain,"beta_uh");
1894  beta_uh_dry            = get_python_double(domain,"beta_uh_dry");
1895  beta_vh                = get_python_double(domain,"beta_vh");
1896  beta_vh_dry            = get_python_double(domain,"beta_vh_dry"); 
1897
1898  minimum_allowed_height = get_python_double(domain,"minimum_allowed_height");
1899  epsilon                = get_python_double(domain,"epsilon");
1900
1901  number_of_elements = stage_centroid_values -> dimensions[0]; 
1902
1903  // Call underlying computational routine
1904  e = _extrapolate_second_order_sw(number_of_elements,
1905                   epsilon,
1906                   minimum_allowed_height,
1907                   beta_w,
1908                   beta_w_dry,
1909                   beta_uh,
1910                   beta_uh_dry,
1911                   beta_vh,
1912                   beta_vh_dry,
1913                   (long*) surrogate_neighbours -> data,
1914                   (long*) number_of_boundaries -> data,
1915                   (double*) centroid_coordinates -> data,
1916                   (double*) stage_centroid_values -> data,
1917                   (double*) xmom_centroid_values -> data,
1918                   (double*) ymom_centroid_values -> data,
1919                   (double*) elevation_centroid_values -> data,
1920                   (double*) vertex_coordinates -> data,
1921                   (double*) stage_vertex_values -> data,
1922                   (double*) xmom_vertex_values -> data,
1923                   (double*) ymom_vertex_values -> data,
1924                   (double*) elevation_vertex_values -> data,
1925                   optimise_dry_cells);
1926  if (e == -1) {
1927    // Use error string set inside computational routine
1928    return NULL;
1929  }                   
1930 
1931 
1932  return Py_BuildValue("");
1933 
1934}// extrapolate_second-order_sw
1935
1936
1937
1938
1939PyObject *rotate(PyObject *self, PyObject *args, PyObject *kwargs) {
1940  //
1941  // r = rotate(q, normal, direction=1)
1942  //
1943  // Where q is assumed to be a Float numeric array of length 3 and
1944  // normal a Float numeric array of length 2.
1945
1946  // FIXME(Ole): I don't think this is used anymore
1947
1948  PyObject *Q, *Normal;
1949  PyArrayObject *q, *r, *normal;
1950
1951  static char *argnames[] = {"q", "normal", "direction", NULL};
1952  int dimensions[1], i, direction=1;
1953  double n1, n2;
1954
1955  // Convert Python arguments to C
1956  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OO|i", argnames,
1957                   &Q, &Normal, &direction)) {
1958    report_python_error(AT, "could not parse input arguments");
1959    return NULL;
1960  } 
1961
1962  // Input checks (convert sequences into numeric arrays)
1963  q = (PyArrayObject *)
1964    PyArray_ContiguousFromObject(Q, PyArray_DOUBLE, 0, 0);
1965  normal = (PyArrayObject *)
1966    PyArray_ContiguousFromObject(Normal, PyArray_DOUBLE, 0, 0);
1967
1968
1969  if (normal -> dimensions[0] != 2) {
1970    report_python_error(AT, "Normal vector must have 2 components");
1971    return NULL;
1972  }
1973
1974  // Allocate space for return vector r (don't DECREF)
1975  dimensions[0] = 3;
1976  r = (PyArrayObject *) anuga_FromDims(1, dimensions, PyArray_DOUBLE);
1977
1978  // Copy
1979  for (i=0; i<3; i++) {
1980    ((double *) (r -> data))[i] = ((double *) (q -> data))[i];
1981  }
1982
1983  // Get normal and direction
1984  n1 = ((double *) normal -> data)[0];
1985  n2 = ((double *) normal -> data)[1];
1986  if (direction == -1) n2 = -n2;
1987
1988  // Rotate
1989  _rotate((double *) r -> data, n1, n2);
1990
1991  // Release numeric arrays
1992  Py_DECREF(q);
1993  Py_DECREF(normal);
1994
1995  // Return result using PyArray to avoid memory leak
1996  return PyArray_Return(r);
1997}
1998
1999
2000// Computational function for flux computation
2001
2002double _compute_fluxes_central(int number_of_elements,
2003        double timestep,
2004        double epsilon,
2005        double H0,
2006        double g,
2007        long* neighbours,
2008        long* neighbour_edges,
2009        double* normals,
2010        double* edgelengths,
2011        double* radii,
2012        double* areas,
2013        long* tri_full_flag,
2014        double* stage_edge_values,
2015        double* xmom_edge_values,
2016        double* ymom_edge_values,
2017        double* bed_edge_values,
2018        double* stage_boundary_values,
2019        double* xmom_boundary_values,
2020        double* ymom_boundary_values,
2021        double* stage_explicit_update,
2022        double* xmom_explicit_update,
2023        double* ymom_explicit_update,
2024        long* already_computed_flux,
2025        double* max_speed_array,
2026        int optimise_dry_cells) {
2027    // Local variables
2028    double max_speed, length, inv_area, zl, zr;
2029    double h0 = H0*H0; // This ensures a good balance when h approaches H0.
2030
2031    double limiting_threshold = 10 * H0; // Avoid applying limiter below this
2032    // threshold for performance reasons.
2033    // See ANUGA manual under flux limiting
2034    int k, i, m, n;
2035    int ki, nm = 0, ki2; // Index shorthands
2036
2037
2038    // Workspace (making them static actually made function slightly slower (Ole))
2039    double ql[3], qr[3], edgeflux[3]; // Work array for summing up fluxes
2040
2041    static long call = 1; // Static local variable flagging already computed flux
2042
2043    // Start computation
2044    call++; // Flag 'id' of flux calculation for this timestep
2045
2046    // Set explicit_update to zero for all conserved_quantities.
2047    // This assumes compute_fluxes called before forcing terms
2048    memset((char*) stage_explicit_update, 0, number_of_elements * sizeof (double));
2049    memset((char*) xmom_explicit_update, 0, number_of_elements * sizeof (double));
2050    memset((char*) ymom_explicit_update, 0, number_of_elements * sizeof (double));
2051
2052    // For all triangles
2053    for (k = 0; k < number_of_elements; k++) {
2054        // Loop through neighbours and compute edge flux for each
2055        for (i = 0; i < 3; i++) {
2056            ki = k * 3 + i; // Linear index to edge i of triangle k
2057
2058            if (already_computed_flux[ki] == call) {
2059                // We've already computed the flux across this edge
2060                continue;
2061            }
2062
2063            // Get left hand side values from triangle k, edge i
2064            ql[0] = stage_edge_values[ki];
2065            ql[1] = xmom_edge_values[ki];
2066            ql[2] = ymom_edge_values[ki];
2067            zl = bed_edge_values[ki];
2068
2069            // Get right hand side values either from neighbouring triangle
2070            // or from boundary array (Quantities at neighbour on nearest face).
2071            n = neighbours[ki];
2072            if (n < 0) {
2073                // Neighbour is a boundary condition
2074                m = -n - 1; // Convert negative flag to boundary index
2075
2076                qr[0] = stage_boundary_values[m];
2077                qr[1] = xmom_boundary_values[m];
2078                qr[2] = ymom_boundary_values[m];
2079                zr = zl; // Extend bed elevation to boundary
2080            }
2081            else {
2082                // Neighbour is a real triangle
2083                m = neighbour_edges[ki];
2084                nm = n * 3 + m; // Linear index (triangle n, edge m)
2085
2086                qr[0] = stage_edge_values[nm];
2087                qr[1] = xmom_edge_values[nm];
2088                qr[2] = ymom_edge_values[nm];
2089                zr = bed_edge_values[nm];
2090            }
2091
2092            // Now we have values for this edge - both from left and right side.
2093
2094            if (optimise_dry_cells) {
2095                // Check if flux calculation is necessary across this edge
2096                // This check will exclude dry cells.
2097                // This will also optimise cases where zl != zr as
2098                // long as both are dry
2099
2100                if (fabs(ql[0] - zl) < epsilon &&
2101                        fabs(qr[0] - zr) < epsilon) {
2102                    // Cell boundary is dry
2103
2104                    already_computed_flux[ki] = call; // #k Done
2105                    if (n >= 0) {
2106                        already_computed_flux[nm] = call; // #n Done
2107                    }
2108
2109                    max_speed = 0.0;
2110                    continue;
2111                }
2112            }
2113
2114
2115            if (fabs(zl-zr)>1.0e-10) {
2116                report_python_error(AT,"Discontinuous Elevation");
2117                return 0.0;
2118            }
2119           
2120            // Outward pointing normal vector (domain.normals[k, 2*i:2*i+2])
2121            ki2 = 2 * ki; //k*6 + i*2
2122
2123            // Edge flux computation (triangle k, edge i)
2124            _flux_function_central(ql, qr, zl, zr,
2125                    normals[ki2], normals[ki2 + 1],
2126                    epsilon, h0, limiting_threshold, g,
2127                    edgeflux, &max_speed);
2128
2129
2130            // Multiply edgeflux by edgelength
2131            length = edgelengths[ki];
2132            edgeflux[0] *= length;
2133            edgeflux[1] *= length;
2134            edgeflux[2] *= length;
2135
2136
2137            // Update triangle k with flux from edge i
2138            stage_explicit_update[k] -= edgeflux[0];
2139            xmom_explicit_update[k] -= edgeflux[1];
2140            ymom_explicit_update[k] -= edgeflux[2];
2141
2142            already_computed_flux[ki] = call; // #k Done
2143
2144
2145            // Update neighbour n with same flux but reversed sign
2146            if (n >= 0) {
2147                stage_explicit_update[n] += edgeflux[0];
2148                xmom_explicit_update[n] += edgeflux[1];
2149                ymom_explicit_update[n] += edgeflux[2];
2150
2151                already_computed_flux[nm] = call; // #n Done
2152            }
2153
2154            // Update timestep based on edge i and possibly neighbour n
2155            if (tri_full_flag[k] == 1) {
2156                if (max_speed > epsilon) {
2157                    // Apply CFL condition for triangles joining this edge (triangle k and triangle n)
2158
2159                    // CFL for triangle k
2160                    timestep = min(timestep, radii[k] / max_speed);
2161
2162                    if (n >= 0) {
2163                        // Apply CFL condition for neigbour n (which is on the ith edge of triangle k)
2164                        timestep = min(timestep, radii[n] / max_speed);
2165                    }
2166
2167                    // Ted Rigby's suggested less conservative version
2168                    //if (n>=0) {
2169                    //  timestep = min(timestep, (radii[k]+radii[n])/max_speed);
2170                    //} else {
2171                    //  timestep = min(timestep, radii[k]/max_speed);
2172                    // }
2173                }
2174            }
2175
2176        } // End edge i (and neighbour n)
2177
2178
2179        // Normalise triangle k by area and store for when all conserved
2180        // quantities get updated
2181        inv_area = 1.0 / areas[k];
2182        stage_explicit_update[k] *= inv_area;
2183        xmom_explicit_update[k] *= inv_area;
2184        ymom_explicit_update[k] *= inv_area;
2185
2186
2187        // Keep track of maximal speeds
2188        max_speed_array[k] = max_speed;
2189
2190    } // End triangle k
2191
2192    return timestep;
2193}
2194
2195//=========================================================================
2196// Python Glue
2197//=========================================================================
2198
2199PyObject *compute_fluxes_ext_central_new(PyObject *self, PyObject *args) {
2200  /*Compute all fluxes and the timestep suitable for all volumes
2201    in domain.
2202
2203    Compute total flux for each conserved quantity using "flux_function_central"
2204
2205    Fluxes across each edge are scaled by edgelengths and summed up
2206    Resulting flux is then scaled by area and stored in
2207    explicit_update for each of the three conserved quantities
2208    stage, xmomentum and ymomentum
2209
2210    The maximal allowable speed computed by the flux_function for each volume
2211    is converted to a timestep that must not be exceeded. The minimum of
2212    those is computed as the next overall timestep.
2213
2214    Python call:
2215    timestep = compute_fluxes(timestep, domain, stage, xmom, ymom, bed)
2216
2217
2218    Post conditions:
2219      domain.explicit_update is reset to computed flux values
2220      returns timestep which is the largest step satisfying all volumes.
2221
2222
2223  */
2224
2225    PyObject
2226    *domain,
2227    *stage, 
2228    *xmom, 
2229    *ymom, 
2230    *bed;
2231
2232    PyArrayObject
2233    *neighbours, 
2234    *neighbour_edges,
2235    *normals, 
2236    *edgelengths, 
2237    *radii, 
2238    *areas,
2239    *tri_full_flag,
2240    *stage_edge_values,
2241    *xmom_edge_values,
2242    *ymom_edge_values,
2243    *bed_edge_values,
2244    *stage_boundary_values,
2245    *xmom_boundary_values,
2246    *ymom_boundary_values,
2247    *stage_explicit_update,
2248    *xmom_explicit_update,
2249    *ymom_explicit_update,
2250    *already_computed_flux, //Tracks whether the flux across an edge has already been computed
2251    *max_speed_array; //Keeps track of max speeds for each triangle
2252
2253   
2254    double timestep, epsilon, H0, g;
2255    int optimise_dry_cells;
2256   
2257    // Convert Python arguments to C
2258    if (!PyArg_ParseTuple(args, "dOOOO", &timestep, &domain, &stage, &xmom, &ymom, &bed )) {
2259      report_python_error(AT, "could not parse input arguments");
2260      return NULL;
2261    }
2262
2263    epsilon           = get_python_double(domain,"epsilon");
2264    H0                = get_python_double(domain,"H0");
2265    g                 = get_python_double(domain,"g");
2266    optimise_dry_cells = get_python_integer(domain,"optimse_dry_cells");
2267   
2268    neighbours        = get_consecutive_array(domain, "neighbours");
2269    neighbour_edges   = get_consecutive_array(domain, "neighbour_edges"); 
2270    normals           = get_consecutive_array(domain, "normals");
2271    edgelengths       = get_consecutive_array(domain, "edge_lengths");   
2272    radii             = get_consecutive_array(domain, "radii");   
2273    areas             = get_consecutive_array(domain, "areas");   
2274    tri_full_flag     = get_consecutive_array(domain, "tri_full_flag");
2275    already_computed_flux  = get_consecutive_array(domain, "already_computed_flux");
2276    max_speed_array   = get_consecutive_array(domain, "max_speed");
2277   
2278    stage_edge_values = get_consecutive_array(stage, "edge_values");   
2279    xmom_edge_values  = get_consecutive_array(xmom, "edge_values");   
2280    ymom_edge_values  = get_consecutive_array(ymom, "edge_values"); 
2281    bed_edge_values   = get_consecutive_array(bed, "edge_values");   
2282
2283    stage_boundary_values = get_consecutive_array(stage, "boundary_values");   
2284    xmom_boundary_values  = get_consecutive_array(xmom, "boundary_values");   
2285    ymom_boundary_values  = get_consecutive_array(ymom, "boundary_values"); 
2286
2287    stage_explicit_update = get_consecutive_array(stage, "explicit_update");   
2288    xmom_explicit_update  = get_consecutive_array(xmom, "explicit_update");   
2289    ymom_explicit_update  = get_consecutive_array(ymom, "explicit_update"); 
2290
2291
2292  int number_of_elements = stage_edge_values -> dimensions[0];
2293
2294  // Call underlying flux computation routine and update
2295  // the explicit update arrays
2296  timestep = _compute_fluxes_central(number_of_elements,
2297                     timestep,
2298                     epsilon,
2299                     H0,
2300                     g,
2301                     (long*) neighbours -> data,
2302                     (long*) neighbour_edges -> data,
2303                     (double*) normals -> data,
2304                     (double*) edgelengths -> data, 
2305                     (double*) radii -> data, 
2306                     (double*) areas -> data,
2307                     (long*) tri_full_flag -> data,
2308                     (double*) stage_edge_values -> data,
2309                     (double*) xmom_edge_values -> data,
2310                     (double*) ymom_edge_values -> data,
2311                     (double*) bed_edge_values -> data,
2312                     (double*) stage_boundary_values -> data,
2313                     (double*) xmom_boundary_values -> data,
2314                     (double*) ymom_boundary_values -> data,
2315                     (double*) stage_explicit_update -> data,
2316                     (double*) xmom_explicit_update -> data,
2317                     (double*) ymom_explicit_update -> data,
2318                     (long*) already_computed_flux -> data,
2319                     (double*) max_speed_array -> data,
2320                     optimise_dry_cells);
2321
2322  Py_DECREF(neighbours);
2323  Py_DECREF(neighbour_edges);
2324  Py_DECREF(normals);
2325  Py_DECREF(edgelengths);
2326  Py_DECREF(radii);
2327  Py_DECREF(areas);
2328  Py_DECREF(tri_full_flag);
2329  Py_DECREF(already_computed_flux);
2330  Py_DECREF(max_speed_array);
2331  Py_DECREF(stage_edge_values);
2332  Py_DECREF(xmom_edge_values);
2333  Py_DECREF(ymom_edge_values);
2334  Py_DECREF(bed_edge_values);
2335  Py_DECREF(stage_boundary_values);
2336  Py_DECREF(xmom_boundary_values);
2337  Py_DECREF(ymom_boundary_values);
2338  Py_DECREF(stage_explicit_update);
2339  Py_DECREF(xmom_explicit_update);
2340  Py_DECREF(ymom_explicit_update);
2341
2342 
2343  // Return updated flux timestep
2344  return Py_BuildValue("d", timestep);
2345}
2346
2347
2348
2349
2350
2351
2352PyObject *compute_fluxes_ext_central(PyObject *self, PyObject *args) {
2353  /*Compute all fluxes and the timestep suitable for all volumes
2354    in domain.
2355
2356    Compute total flux for each conserved quantity using "flux_function_central"
2357
2358    Fluxes across each edge are scaled by edgelengths and summed up
2359    Resulting flux is then scaled by area and stored in
2360    explicit_update for each of the three conserved quantities
2361    stage, xmomentum and ymomentum
2362
2363    The maximal allowable speed computed by the flux_function for each volume
2364    is converted to a timestep that must not be exceeded. The minimum of
2365    those is computed as the next overall timestep.
2366
2367    Python call:
2368    domain.timestep = compute_fluxes(timestep,
2369                                     domain.epsilon,
2370                     domain.H0,
2371                                     domain.g,
2372                                     domain.neighbours,
2373                                     domain.neighbour_edges,
2374                                     domain.normals,
2375                                     domain.edgelengths,
2376                                     domain.radii,
2377                                     domain.areas,
2378                                     tri_full_flag,
2379                                     Stage.edge_values,
2380                                     Xmom.edge_values,
2381                                     Ymom.edge_values,
2382                                     Bed.edge_values,
2383                                     Stage.boundary_values,
2384                                     Xmom.boundary_values,
2385                                     Ymom.boundary_values,
2386                                     Stage.explicit_update,
2387                                     Xmom.explicit_update,
2388                                     Ymom.explicit_update,
2389                                     already_computed_flux,
2390                     optimise_dry_cells)                     
2391
2392
2393    Post conditions:
2394      domain.explicit_update is reset to computed flux values
2395      domain.timestep is set to the largest step satisfying all volumes.
2396
2397
2398  */
2399
2400
2401  PyArrayObject *neighbours, *neighbour_edges,
2402    *normals, *edgelengths, *radii, *areas,
2403    *tri_full_flag,
2404    *stage_edge_values,
2405    *xmom_edge_values,
2406    *ymom_edge_values,
2407    *bed_edge_values,
2408    *stage_boundary_values,
2409    *xmom_boundary_values,
2410    *ymom_boundary_values,
2411    *stage_explicit_update,
2412    *xmom_explicit_update,
2413    *ymom_explicit_update,
2414    *already_computed_flux, //Tracks whether the flux across an edge has already been computed
2415    *max_speed_array; //Keeps track of max speeds for each triangle
2416
2417   
2418  double timestep, epsilon, H0, g;
2419  int optimise_dry_cells;
2420   
2421  // Convert Python arguments to C
2422  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOi",
2423            &timestep,
2424            &epsilon,
2425            &H0,
2426            &g,
2427            &neighbours,
2428            &neighbour_edges,
2429            &normals,
2430            &edgelengths, &radii, &areas,
2431            &tri_full_flag,
2432            &stage_edge_values,
2433            &xmom_edge_values,
2434            &ymom_edge_values,
2435            &bed_edge_values,
2436            &stage_boundary_values,
2437            &xmom_boundary_values,
2438            &ymom_boundary_values,
2439            &stage_explicit_update,
2440            &xmom_explicit_update,
2441            &ymom_explicit_update,
2442            &already_computed_flux,
2443            &max_speed_array,
2444            &optimise_dry_cells)) {
2445    report_python_error(AT, "could not parse input arguments");
2446    return NULL;
2447  }
2448
2449  // check that numpy array objects arrays are C contiguous memory
2450  CHECK_C_CONTIG(neighbours);
2451  CHECK_C_CONTIG(neighbour_edges);
2452  CHECK_C_CONTIG(normals);
2453  CHECK_C_CONTIG(edgelengths);
2454  CHECK_C_CONTIG(radii);
2455  CHECK_C_CONTIG(areas);
2456  CHECK_C_CONTIG(tri_full_flag);
2457  CHECK_C_CONTIG(stage_edge_values);
2458  CHECK_C_CONTIG(xmom_edge_values);
2459  CHECK_C_CONTIG(ymom_edge_values);
2460  CHECK_C_CONTIG(bed_edge_values);
2461  CHECK_C_CONTIG(stage_boundary_values);
2462  CHECK_C_CONTIG(xmom_boundary_values);
2463  CHECK_C_CONTIG(ymom_boundary_values);
2464  CHECK_C_CONTIG(stage_explicit_update);
2465  CHECK_C_CONTIG(xmom_explicit_update);
2466  CHECK_C_CONTIG(ymom_explicit_update);
2467  CHECK_C_CONTIG(already_computed_flux);
2468  CHECK_C_CONTIG(max_speed_array);
2469
2470  int number_of_elements = stage_edge_values -> dimensions[0];
2471
2472  // Call underlying flux computation routine and update
2473  // the explicit update arrays
2474  timestep = _compute_fluxes_central(number_of_elements,
2475                     timestep,
2476                     epsilon,
2477                     H0,
2478                     g,
2479                     (long*) neighbours -> data,
2480                     (long*) neighbour_edges -> data,
2481                     (double*) normals -> data,
2482                     (double*) edgelengths -> data, 
2483                     (double*) radii -> data, 
2484                     (double*) areas -> data,
2485                     (long*) tri_full_flag -> data,
2486                     (double*) stage_edge_values -> data,
2487                     (double*) xmom_edge_values -> data,
2488                     (double*) ymom_edge_values -> data,
2489                     (double*) bed_edge_values -> data,
2490                     (double*) stage_boundary_values -> data,
2491                     (double*) xmom_boundary_values -> data,
2492                     (double*) ymom_boundary_values -> data,
2493                     (double*) stage_explicit_update -> data,
2494                     (double*) xmom_explicit_update -> data,
2495                     (double*) ymom_explicit_update -> data,
2496                     (long*) already_computed_flux -> data,
2497                     (double*) max_speed_array -> data,
2498                     optimise_dry_cells);
2499 
2500  // Return updated flux timestep
2501  return Py_BuildValue("d", timestep);
2502}
2503
2504
2505
2506
2507
2508PyObject *compute_fluxes_ext_kinetic(PyObject *self, PyObject *args) {
2509  /*Compute all fluxes and the timestep suitable for all volumes
2510    in domain.
2511
2512    THIS IS AN EXPERIMENTAL FUNCTION - NORMALLY flux_function_central IS USED.
2513   
2514    Compute total flux for each conserved quantity using "flux_function_kinetic"
2515
2516    Fluxes across each edge are scaled by edgelengths and summed up
2517    Resulting flux is then scaled by area and stored in
2518    explicit_update for each of the three conserved quantities
2519    stage, xmomentum and ymomentum
2520
2521    The maximal allowable speed computed by the flux_function for each volume
2522    is converted to a timestep that must not be exceeded. The minimum of
2523    those is computed as the next overall timestep.
2524
2525    Python call:
2526    domain.timestep = compute_fluxes(timestep,
2527                                     domain.epsilon,
2528                                     domain.H0,
2529                                     domain.g,
2530                                     domain.neighbours,
2531                                     domain.neighbour_edges,
2532                                     domain.normals,
2533                                     domain.edgelengths,
2534                                     domain.radii,
2535                                     domain.areas,
2536                                     Stage.edge_values,
2537                                     Xmom.edge_values,
2538                                     Ymom.edge_values,
2539                                     Bed.edge_values,
2540                                     Stage.boundary_values,
2541                                     Xmom.boundary_values,
2542                                     Ymom.boundary_values,
2543                                     Stage.explicit_update,
2544                                     Xmom.explicit_update,
2545                                     Ymom.explicit_update,
2546                                     already_computed_flux)
2547
2548
2549    Post conditions:
2550      domain.explicit_update is reset to computed flux values
2551      domain.timestep is set to the largest step satisfying all volumes.
2552
2553
2554  */
2555
2556
2557  PyArrayObject *neighbours, *neighbour_edges,
2558    *normals, *edgelengths, *radii, *areas,
2559    *tri_full_flag,
2560    *stage_edge_values,
2561    *xmom_edge_values,
2562    *ymom_edge_values,
2563    *bed_edge_values,
2564    *stage_boundary_values,
2565    *xmom_boundary_values,
2566    *ymom_boundary_values,
2567    *stage_explicit_update,
2568    *xmom_explicit_update,
2569    *ymom_explicit_update,
2570    *already_computed_flux; // Tracks whether the flux across an edge has already been computed
2571
2572
2573  // Local variables
2574  double timestep, max_speed, epsilon, g, H0;
2575  double normal[2], ql[3], qr[3], zl, zr;
2576  double edgeflux[3]; //Work arrays for summing up fluxes
2577
2578  int number_of_elements, k, i, m, n;
2579  int ki, nm=0, ki2; //Index shorthands
2580  static long call=1;
2581
2582
2583  // Convert Python arguments to C
2584  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOO",
2585            &timestep,
2586            &epsilon,
2587            &H0,
2588            &g,
2589            &neighbours,
2590            &neighbour_edges,
2591            &normals,
2592            &edgelengths, &radii, &areas,
2593            &tri_full_flag,
2594            &stage_edge_values,
2595            &xmom_edge_values,
2596            &ymom_edge_values,
2597            &bed_edge_values,
2598            &stage_boundary_values,
2599            &xmom_boundary_values,
2600            &ymom_boundary_values,
2601            &stage_explicit_update,
2602            &xmom_explicit_update,
2603            &ymom_explicit_update,
2604            &already_computed_flux)) {
2605    report_python_error(AT, "could not parse input arguments");
2606    return NULL;
2607  }
2608  number_of_elements = stage_edge_values -> dimensions[0];
2609  call++;//a static local variable to which already_computed_flux is compared
2610  //set explicit_update to zero for all conserved_quantities.
2611  //This assumes compute_fluxes called before forcing terms
2612  for (k=0; k<number_of_elements; k++) {
2613    ((double *) stage_explicit_update -> data)[k]=0.0;
2614    ((double *) xmom_explicit_update -> data)[k]=0.0;
2615    ((double *) ymom_explicit_update -> data)[k]=0.0;
2616  }
2617  //Loop through neighbours and compute edge flux for each
2618  for (k=0; k<number_of_elements; k++) {
2619    for (i=0; i<3; i++) {
2620      ki = k*3+i;
2621      if (((long *) already_computed_flux->data)[ki]==call)//we've already computed the flux across this edge
2622    continue;
2623      ql[0] = ((double *) stage_edge_values -> data)[ki];
2624      ql[1] = ((double *) xmom_edge_values -> data)[ki];
2625      ql[2] = ((double *) ymom_edge_values -> data)[ki];
2626      zl =    ((double *) bed_edge_values -> data)[ki];
2627
2628      //Quantities at neighbour on nearest face
2629      n = ((long *) neighbours -> data)[ki];
2630      if (n < 0) {
2631    m = -n-1; //Convert negative flag to index
2632    qr[0] = ((double *) stage_boundary_values -> data)[m];
2633    qr[1] = ((double *) xmom_boundary_values -> data)[m];
2634    qr[2] = ((double *) ymom_boundary_values -> data)[m];
2635    zr = zl; //Extend bed elevation to boundary
2636      } else {
2637    m = ((long *) neighbour_edges -> data)[ki];
2638    nm = n*3+m;
2639    qr[0] = ((double *) stage_edge_values -> data)[nm];
2640    qr[1] = ((double *) xmom_edge_values -> data)[nm];
2641    qr[2] = ((double *) ymom_edge_values -> data)[nm];
2642    zr =    ((double *) bed_edge_values -> data)[nm];
2643      }
2644      // Outward pointing normal vector
2645      // normal = domain.normals[k, 2*i:2*i+2]
2646      ki2 = 2*ki; //k*6 + i*2
2647      normal[0] = ((double *) normals -> data)[ki2];
2648      normal[1] = ((double *) normals -> data)[ki2+1];
2649      //Edge flux computation
2650      flux_function_kinetic(ql, qr, zl, zr,
2651            normal[0], normal[1],
2652            epsilon, H0, g,
2653            edgeflux, &max_speed);
2654      //update triangle k
2655      ((long *) already_computed_flux->data)[ki]=call;
2656      ((double *) stage_explicit_update -> data)[k] -= edgeflux[0]*((double *) edgelengths -> data)[ki];
2657      ((double *) xmom_explicit_update -> data)[k] -= edgeflux[1]*((double *) edgelengths -> data)[ki];
2658      ((double *) ymom_explicit_update -> data)[k] -= edgeflux[2]*((double *) edgelengths -> data)[ki];
2659      //update the neighbour n
2660      if (n>=0){
2661    ((long *) already_computed_flux->data)[nm]=call;
2662    ((double *) stage_explicit_update -> data)[n] += edgeflux[0]*((double *) edgelengths -> data)[nm];
2663    ((double *) xmom_explicit_update -> data)[n] += edgeflux[1]*((double *) edgelengths -> data)[nm];
2664    ((double *) ymom_explicit_update -> data)[n] += edgeflux[2]*((double *) edgelengths -> data)[nm];
2665      }
2666      ///for (j=0; j<3; j++) {
2667    ///flux[j] -= edgeflux[j]*((double *) edgelengths -> data)[ki];
2668    ///}
2669    //Update timestep
2670    //timestep = min(timestep, domain.radii[k]/max_speed)
2671    //FIXME: SR Add parameter for CFL condition
2672    if ( ((long *) tri_full_flag -> data)[k] == 1) {
2673        if (max_speed > epsilon) {
2674            timestep = min(timestep, ((double *) radii -> data)[k]/max_speed);
2675            //maxspeed in flux_function is calculated as max(|u+a|,|u-a|)
2676            if (n>=0)
2677                timestep = min(timestep, ((double *) radii -> data)[n]/max_speed);
2678        }
2679    }
2680    } // end for i
2681    //Normalise by area and store for when all conserved
2682    //quantities get updated
2683    ((double *) stage_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
2684    ((double *) xmom_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
2685    ((double *) ymom_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
2686  } //end for k
2687  return Py_BuildValue("d", timestep);
2688}
2689
2690PyObject *protect(PyObject *self, PyObject *args) {
2691  //
2692  //    protect(minimum_allowed_height, maximum_allowed_speed, wc, zc, xmomc, ymomc)
2693
2694
2695  PyArrayObject
2696  *wc,            //Stage at centroids
2697  *zc,            //Elevation at centroids
2698  *xmomc,         //Momentums at centroids
2699  *ymomc;
2700
2701
2702  int N;
2703  double minimum_allowed_height, maximum_allowed_speed, epsilon;
2704
2705  // Convert Python arguments to C
2706  if (!PyArg_ParseTuple(args, "dddOOOO",
2707            &minimum_allowed_height,
2708            &maximum_allowed_speed,
2709            &epsilon,
2710            &wc, &zc, &xmomc, &ymomc)) {
2711    report_python_error(AT, "could not parse input arguments");
2712    return NULL;
2713  } 
2714
2715  N = wc -> dimensions[0];
2716
2717  _protect(N,
2718       minimum_allowed_height,
2719       maximum_allowed_speed,
2720       epsilon,
2721       (double*) wc -> data,
2722       (double*) zc -> data,
2723       (double*) xmomc -> data,
2724       (double*) ymomc -> data);
2725
2726  return Py_BuildValue("");
2727}
2728
2729
2730
2731PyObject *balance_deep_and_shallow(PyObject *self, PyObject *args) {
2732  // Compute linear combination between stage as computed by
2733  // gradient-limiters limiting using w, and stage computed by
2734  // gradient-limiters limiting using h (h-limiter).
2735  // The former takes precedence when heights are large compared to the
2736  // bed slope while the latter takes precedence when heights are
2737  // relatively small.  Anything in between is computed as a balanced
2738  // linear combination in order to avoid numerical disturbances which
2739  // would otherwise appear as a result of hard switching between
2740  // modes.
2741  //
2742  //    balance_deep_and_shallow(wc, zc, wv, zv,
2743  //                             xmomc, ymomc, xmomv, ymomv)
2744
2745
2746  PyArrayObject
2747    *wc,            //Stage at centroids
2748    *zc,            //Elevation at centroids
2749    *wv,            //Stage at vertices
2750    *zv,            //Elevation at vertices
2751    *hvbar,         //h-Limited depths at vertices
2752    *xmomc,         //Momentums at centroids and vertices
2753    *ymomc,
2754    *xmomv,
2755    *ymomv;
2756 
2757  PyObject *domain, *Tmp;
2758   
2759  double alpha_balance = 2.0;
2760  double H0;
2761
2762  int N, tight_slope_limiters, use_centroid_velocities; //, err;
2763
2764  // Convert Python arguments to C
2765  if (!PyArg_ParseTuple(args, "OOOOOOOOOO",
2766            &domain,
2767            &wc, &zc, 
2768            &wv, &zv, &hvbar,
2769            &xmomc, &ymomc, &xmomv, &ymomv)) {
2770    report_python_error(AT, "could not parse input arguments");
2771    return NULL;
2772  } 
2773     
2774     
2775  // FIXME (Ole): I tested this without GetAttrString and got time down
2776  // marginally from 4.0s to 3.8s. Consider passing everything in
2777  // through ParseTuple and profile.
2778 
2779  // Pull out parameters
2780  Tmp = PyObject_GetAttrString(domain, "alpha_balance");
2781  if (!Tmp) {
2782    report_python_error(AT, "could not obtain object alpha_balance from domain");
2783    return NULL;
2784  } 
2785  alpha_balance = PyFloat_AsDouble(Tmp);
2786  Py_DECREF(Tmp);
2787
2788 
2789  Tmp = PyObject_GetAttrString(domain, "H0");
2790  if (!Tmp) {
2791    report_python_error(AT, "could not obtain object H0 from domain");
2792    return NULL;
2793  } 
2794  H0 = PyFloat_AsDouble(Tmp);
2795  Py_DECREF(Tmp);
2796
2797 
2798  Tmp = PyObject_GetAttrString(domain, "tight_slope_limiters");
2799  if (!Tmp) {
2800    report_python_error(AT, "could not obtain object tight_slope_limiters from domain");
2801    return NULL;
2802  } 
2803  tight_slope_limiters = PyInt_AsLong(Tmp);
2804  Py_DECREF(Tmp);
2805 
2806 
2807  Tmp = PyObject_GetAttrString(domain, "use_centroid_velocities");
2808  if (!Tmp) {
2809    report_python_error(AT, "could not obtain object use_centroid_velocities from domain");
2810    return NULL;
2811  } 
2812  use_centroid_velocities = PyInt_AsLong(Tmp);
2813  Py_DECREF(Tmp);
2814 
2815 
2816 
2817  N = wc -> dimensions[0];
2818  _balance_deep_and_shallow(N,
2819                (double*) wc -> data,
2820                (double*) zc -> data,
2821                (double*) wv -> data,
2822                (double*) zv -> data,
2823                (double*) hvbar -> data,
2824                (double*) xmomc -> data,
2825                (double*) ymomc -> data,
2826                (double*) xmomv -> data,
2827                (double*) ymomv -> data,
2828                H0,
2829                (int) tight_slope_limiters,
2830                (int) use_centroid_velocities,             
2831                alpha_balance);
2832
2833
2834  return Py_BuildValue("");
2835}
2836
2837
2838
2839
2840PyObject *assign_windfield_values(PyObject *self, PyObject *args) {
2841  //
2842  //      assign_windfield_values(xmom_update, ymom_update,
2843  //                              s_vec, phi_vec, self.const)
2844
2845
2846
2847  PyArrayObject   //(one element per triangle)
2848  *s_vec,         //Speeds
2849  *phi_vec,       //Bearings
2850  *xmom_update,   //Momentum updates
2851  *ymom_update;
2852
2853
2854  int N;
2855  double cw;
2856
2857  // Convert Python arguments to C
2858  if (!PyArg_ParseTuple(args, "OOOOd",
2859            &xmom_update,
2860            &ymom_update,
2861            &s_vec, &phi_vec,
2862            &cw)) {
2863    report_python_error(AT, "could not parse input arguments");
2864    return NULL;
2865  } 
2866           
2867
2868  N = xmom_update -> dimensions[0];
2869
2870  _assign_wind_field_values(N,
2871       (double*) xmom_update -> data,
2872       (double*) ymom_update -> data,
2873       (double*) s_vec -> data,
2874       (double*) phi_vec -> data,
2875       cw);
2876
2877  return Py_BuildValue("");
2878}
2879
2880
2881
2882
2883//-------------------------------
2884// Method table for python module
2885//-------------------------------
2886static struct PyMethodDef MethodTable[] = {
2887  /* The cast of the function is necessary since PyCFunction values
2888   * only take two PyObject* parameters, and rotate() takes
2889   * three.
2890   */
2891
2892  {"rotate", (PyCFunction)rotate, METH_VARARGS | METH_KEYWORDS, "Print out"},
2893  {"extrapolate_second_order_sw", extrapolate_second_order_sw, METH_VARARGS, "Print out"},
2894  {"compute_fluxes_ext_central", compute_fluxes_ext_central, METH_VARARGS, "Print out"},
2895  {"compute_fluxes_ext_central_new", compute_fluxes_ext_central_new, METH_VARARGS, "Print out"},
2896  {"compute_fluxes_ext_kinetic", compute_fluxes_ext_kinetic, METH_VARARGS, "Print out"},
2897  {"gravity", gravity, METH_VARARGS, "Print out"},
2898  {"manning_friction_old", manning_friction_old, METH_VARARGS, "Print out"},
2899  {"manning_friction_new", manning_friction_new, METH_VARARGS, "Print out"},
2900  {"chezy_friction", chezy_friction, METH_VARARGS, "Print out"},
2901  {"flux_function_central", flux_function_central, METH_VARARGS, "Print out"}, 
2902  {"balance_deep_and_shallow", balance_deep_and_shallow,
2903   METH_VARARGS, "Print out"},
2904  {"protect", protect, METH_VARARGS | METH_KEYWORDS, "Print out"},
2905  {"assign_windfield_values", assign_windfield_values,
2906   METH_VARARGS | METH_KEYWORDS, "Print out"},
2907  {NULL, NULL}
2908};
2909
2910// Module initialisation
2911void initshallow_water_ext(void){
2912  Py_InitModule("shallow_water_ext", MethodTable);
2913
2914  import_array(); // Necessary for handling of NumPY structures
2915}
Note: See TracBrowser for help on using the repository browser.