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

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

Slight conflict with init.py due to addition of hole_tags

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