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

Last change on this file since 7844 was 7585, checked in by ole, 14 years ago

Added better error message in the rare case of negative triangle area

File size: 86.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/*
599void _manning_friction_explicit(double g, double eps, int N,
600               double* w, double* z,
601               double* uh, double* vh,
602               double* eta, double* xmom, double* ymom) {
603
604  int k;
605  double S, h;
606
607  for (k=0; k<N; k++) {
608    if (eta[k] > eps) {
609      h = w[k]-z[k];
610      if (h >= eps) {
611    S = -g * eta[k]*eta[k] * sqrt((uh[k]*uh[k] + vh[k]*vh[k]));
612    S /= pow(h, 7.0/3);      //Expensive (on Ole's home computer)
613    //S /= exp(7.0/3.0*log(h));      //seems to save about 15% over manning_friction
614    //S /= h*h*(1 + h/3.0 - h*h/9.0); //FIXME: Could use a Taylor expansion
615
616
617    //Update momentum
618    xmom[k] += S*uh[k];
619    ymom[k] += S*vh[k];
620      }
621    }
622  }
623}
624*/
625
626
627
628
629double velocity_balance(double uh_i, 
630            double uh,
631            double h_i, 
632            double h, 
633            double alpha,
634            double epsilon) {
635  // Find alpha such that speed at the vertex is within one
636  // order of magnitude of the centroid speed   
637
638  // FIXME(Ole): Work in progress
639 
640  double a, b, estimate;
641  double m=10; // One order of magnitude - allow for velocity deviations at vertices
642
643 
644  printf("alpha = %f, uh_i=%f, uh=%f, h_i=%f, h=%f\n",
645     alpha, uh_i, uh, h_i, h);
646     
647 
648 
649   
650  // Shorthands and determine inequality
651  if (fabs(uh) < epsilon) {
652    a = 1.0e10; // Limit
653  } else {
654    a = fabs(uh_i - uh)/fabs(uh);
655  }
656     
657  if (h < epsilon) {
658    b = 1.0e10; // Limit
659  } else {   
660    b = m*fabs(h_i - h)/h;
661  }
662
663  printf("a %f, b %f\n", a, b); 
664
665  if (a > b) {
666    estimate = (m-1)/(a-b);     
667   
668    printf("Alpha %f, estimate %f\n", 
669       alpha, estimate);   
670       
671    if (alpha < estimate) {
672      printf("Adjusting alpha from %f to %f\n", 
673         alpha, estimate);
674      alpha = estimate;
675    }
676  } else {
677 
678    if (h < h_i) {
679      estimate = (m-1)/(a-b);             
680   
681      printf("Alpha %f, estimate %f\n", 
682         alpha, estimate);   
683       
684      if (alpha < estimate) {
685    printf("Adjusting alpha from %f to %f\n", 
686           alpha, estimate);
687    alpha = estimate;
688      }   
689    }
690    // Always fulfilled as alpha and m-1 are non negative
691  }
692 
693 
694  return alpha;
695}
696
697
698int _balance_deep_and_shallow(int N,
699                  double* wc,
700                  double* zc,
701                  double* wv,
702                  double* zv,
703                  double* hvbar, // Retire this
704                  double* xmomc,
705                  double* ymomc,
706                  double* xmomv,
707                  double* ymomv,
708                  double H0,
709                  int tight_slope_limiters,
710                  int use_centroid_velocities,
711                  double alpha_balance) {
712
713  int k, k3, i; 
714
715  double dz, hmin, alpha, h_diff, hc_k;
716  double epsilon = 1.0e-6; // FIXME: Temporary measure
717  double hv[3]; // Depths at vertices
718  double uc, vc; // Centroid speeds
719
720  // Compute linear combination between w-limited stages and
721  // h-limited stages close to the bed elevation.
722
723  for (k=0; k<N; k++) {
724    // Compute maximal variation in bed elevation
725    // This quantitiy is
726    //     dz = max_i abs(z_i - z_c)
727    // and it is independent of dimension
728    // In the 1d case zc = (z0+z1)/2
729    // In the 2d case zc = (z0+z1+z2)/3
730
731    k3 = 3*k;
732    hc_k = wc[k] - zc[k]; // Centroid value at triangle k
733   
734    dz = 0.0;
735    if (tight_slope_limiters == 0) {     
736      // FIXME: Try with this one precomputed
737      for (i=0; i<3; i++) {
738        dz = max(dz, fabs(zv[k3+i]-zc[k]));
739      }
740    }
741
742    // Calculate depth at vertices (possibly negative here!)
743    hv[0] = wv[k3] -   zv[k3];
744    hv[1] = wv[k3+1] - zv[k3+1];
745    hv[2] = wv[k3+2] - zv[k3+2];       
746   
747    // Calculate minimal depth across all three vertices
748    hmin = min(hv[0], min(hv[1], hv[2]));
749
750    //if (hmin < 0.0 ) {
751    //  printf("hmin = %f\n",hmin);
752    //}
753   
754
755    // Create alpha in [0,1], where alpha==0 means using the h-limited
756    // stage and alpha==1 means using the w-limited stage as
757    // computed by the gradient limiter (both 1st or 2nd order)
758    if (tight_slope_limiters == 0) {     
759      // If hmin > dz/alpha_balance then alpha = 1 and the bed will have no
760      // effect
761      // If hmin < 0 then alpha = 0 reverting to constant height above bed.
762      // The parameter alpha_balance==2 by default
763
764     
765      if (dz > 0.0) {
766        alpha = max( min( alpha_balance*hmin/dz, 1.0), 0.0 );     
767      } else {
768        alpha = 1.0;  // Flat bed
769      }
770      //printf("Using old style limiter\n");
771     
772    } else {
773
774      // Tight Slope Limiter (2007)
775   
776      // Make alpha as large as possible but still ensure that
777      // final depth is positive and that velocities at vertices
778      // are controlled
779   
780      if (hmin < H0) {
781        alpha = 1.0;
782        for (i=0; i<3; i++) {
783         
784          h_diff = hc_k - hv[i];     
785          if (h_diff <= 0) {
786            // Deep water triangle is further away from bed than
787            // shallow water (hbar < h). Any alpha will do
788     
789          } else { 
790            // Denominator is positive which means that we need some of the
791            // h-limited stage.
792           
793            alpha = min(alpha, (hc_k - H0)/h_diff);     
794          }
795        }
796
797        // Ensure alpha in [0,1]
798        if (alpha>1.0) alpha=1.0;
799        if (alpha<0.0) alpha=0.0;
800       
801      } else {
802        // Use w-limited stage exclusively in deeper water.
803        alpha = 1.0;       
804      }
805    }
806   
807   
808    //  Let
809    //
810    //    wvi be the w-limited stage (wvi = zvi + hvi)
811    //    wvi- be the h-limited state (wvi- = zvi + hvi-)
812    //
813    //
814    //  where i=0,1,2 denotes the vertex ids
815    //
816    //  Weighted balance between w-limited and h-limited stage is
817    //
818    //    wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi)
819    //
820    //  It follows that the updated wvi is
821    //    wvi := zvi + (1-alpha)*hvi- + alpha*hvi
822    //
823    //   Momentum is balanced between constant and limited
824
825   
826    if (alpha < 1) {     
827      for (i=0; i<3; i++) {
828       
829        wv[k3+i] = zv[k3+i] + (1-alpha)*hc_k + alpha*hv[i]; 
830
831        // Update momentum at vertices
832        if (use_centroid_velocities == 1) {
833          // This is a simple, efficient and robust option
834          // It uses first order approximation of velocities, but retains
835          // the order used by stage.
836   
837          // Speeds at centroids
838          if (hc_k > epsilon) {
839            uc = xmomc[k]/hc_k;
840            vc = ymomc[k]/hc_k;
841          } else {
842            uc = 0.0;
843            vc = 0.0;
844          }
845         
846          // Vertex momenta guaranteed to be consistent with depth guaranteeing
847          // controlled speed
848          hv[i] = wv[k3+i] - zv[k3+i]; // Recompute (balanced) vertex depth
849          xmomv[k3+i] = uc*hv[i];
850          ymomv[k3+i] = vc*hv[i];
851     
852        } else {
853          // Update momentum as a linear combination of
854          // xmomc and ymomc (shallow) and momentum
855          // from extrapolator xmomv and ymomv (deep).
856          // This assumes that values from xmomv and ymomv have
857          // been established e.g. by the gradient limiter.
858         
859          // FIXME (Ole): I think this should be used with vertex momenta
860          // computed above using centroid_velocities instead of xmomc
861          // and ymomc as they'll be more representative first order
862          // values.
863         
864          xmomv[k3+i] = (1-alpha)*xmomc[k] + alpha*xmomv[k3+i];
865          ymomv[k3+i] = (1-alpha)*ymomc[k] + alpha*ymomv[k3+i];
866         
867        }
868      }
869    }
870  }
871  return 0;
872}
873 
874
875
876
877int _protect(int N,
878         double minimum_allowed_height,
879         double maximum_allowed_speed,
880         double epsilon,
881         double* wc,
882         double* zc,
883         double* xmomc,
884         double* ymomc) {
885
886  int k;
887  double hc;
888  double u, v, reduced_speed;
889
890
891  // Protect against initesimal and negative heights 
892  if (maximum_allowed_speed < epsilon) {
893    for (k=0; k<N; k++) {
894      hc = wc[k] - zc[k];
895
896      if (hc < minimum_allowed_height) {
897       
898        // Set momentum to zero and ensure h is non negative
899        xmomc[k] = 0.0;
900        ymomc[k] = 0.0;
901        if (hc <= 0.0) wc[k] = zc[k];
902      }
903    }
904  } else {
905   
906    // Protect against initesimal and negative heights
907    for (k=0; k<N; k++) {
908      hc = wc[k] - zc[k];
909     
910      if (hc < minimum_allowed_height) {
911
912        //New code: Adjust momentum to guarantee speeds are physical
913        //          ensure h is non negative
914             
915        if (hc <= 0.0) {
916            wc[k] = zc[k];
917        xmomc[k] = 0.0;
918        ymomc[k] = 0.0;
919        } else {
920          //Reduce excessive speeds derived from division by small hc
921          //FIXME (Ole): This may be unnecessary with new slope limiters
922          //in effect.
923         
924          u = xmomc[k]/hc;
925          if (fabs(u) > maximum_allowed_speed) {
926            reduced_speed = maximum_allowed_speed * u/fabs(u);
927            //printf("Speed (u) has been reduced from %.3f to %.3f\n",
928            //   u, reduced_speed);
929            xmomc[k] = reduced_speed * hc;
930          }
931         
932          v = ymomc[k]/hc;
933          if (fabs(v) > maximum_allowed_speed) {
934            reduced_speed = maximum_allowed_speed * v/fabs(v);
935            //printf("Speed (v) has been reduced from %.3f to %.3f\n",
936            //   v, reduced_speed);
937            ymomc[k] = reduced_speed * hc;
938          }
939        }
940      }
941    }
942  }
943  return 0;
944}
945
946
947
948int _assign_wind_field_values(int N,
949                  double* xmom_update,
950                  double* ymom_update,
951                  double* s_vec,
952                  double* phi_vec,
953                  double cw) {
954
955  // Assign windfield values to momentum updates
956
957  int k;
958  double S, s, phi, u, v;
959
960  for (k=0; k<N; k++) {
961
962    s = s_vec[k];
963    phi = phi_vec[k];
964
965    //Convert to radians
966    phi = phi*pi/180;
967
968    //Compute velocity vector (u, v)
969    u = s*cos(phi);
970    v = s*sin(phi);
971
972    //Compute wind stress
973    S = cw * sqrt(u*u + v*v);
974    xmom_update[k] += S*u;
975    ymom_update[k] += S*v;
976  }
977  return 0;
978}
979
980
981
982///////////////////////////////////////////////////////////////////
983// Gateways to Python
984
985
986
987PyObject *flux_function_central(PyObject *self, PyObject *args) {
988  //
989  // Gateway to innermost flux function.
990  // This is only used by the unit tests as the c implementation is
991  // normally called by compute_fluxes in this module.
992
993
994  PyArrayObject *normal, *ql, *qr,  *edgeflux;
995  double g, epsilon, max_speed, H0, zl, zr;
996  double h0, limiting_threshold;
997
998  if (!PyArg_ParseTuple(args, "OOOddOddd",
999            &normal, &ql, &qr, &zl, &zr, &edgeflux,
1000            &epsilon, &g, &H0)) {
1001    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: flux_function_central could not parse input arguments");
1002    return NULL;
1003  }
1004
1005 
1006  h0 = H0*H0; // This ensures a good balance when h approaches H0.
1007              // But evidence suggests that h0 can be as little as
1008              // epsilon!
1009             
1010  limiting_threshold = 10*H0; // Avoid applying limiter below this
1011                              // threshold for performance reasons.
1012                              // See ANUGA manual under flux limiting 
1013 
1014  _flux_function_central((double*) ql -> data, 
1015                         (double*) qr -> data, 
1016                         zl, 
1017                         zr,                         
1018                         ((double*) normal -> data)[0],
1019                         ((double*) normal -> data)[1],         
1020                         epsilon, h0, limiting_threshold,
1021                         g,
1022                         (double*) edgeflux -> data, 
1023                         &max_speed);
1024 
1025  return Py_BuildValue("d", max_speed); 
1026}
1027
1028
1029
1030PyObject *gravity(PyObject *self, PyObject *args) {
1031  //
1032  //  gravity(g, h, v, x, xmom, ymom)
1033  //
1034
1035
1036  PyArrayObject *h, *v, *x, *xmom, *ymom;
1037  int k, N, k3, k6;
1038  double g, avg_h, zx, zy;
1039  double x0, y0, x1, y1, x2, y2, z0, z1, z2;
1040  //double epsilon;
1041
1042  if (!PyArg_ParseTuple(args, "dOOOOO",
1043            &g, &h, &v, &x,
1044            &xmom, &ymom)) {
1045    //&epsilon)) {
1046    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: gravity could not parse input arguments");
1047    return NULL;
1048  }
1049
1050  // check that numpy array objects arrays are C contiguous memory
1051  CHECK_C_CONTIG(h);
1052  CHECK_C_CONTIG(v);
1053  CHECK_C_CONTIG(x);
1054  CHECK_C_CONTIG(xmom);
1055  CHECK_C_CONTIG(ymom);
1056
1057  N = h -> dimensions[0];
1058  for (k=0; k<N; k++) {
1059    k3 = 3*k;  // base index
1060
1061    // Get bathymetry
1062    z0 = ((double*) v -> data)[k3 + 0];
1063    z1 = ((double*) v -> data)[k3 + 1];
1064    z2 = ((double*) v -> data)[k3 + 2];
1065
1066    // Optimise for flat bed
1067    // Note (Ole): This didn't produce measurable speed up.
1068    // Revisit later
1069    //if (fabs(z0-z1)<epsilon && fabs(z1-z2)<epsilon) {
1070    //  continue;
1071    //}
1072
1073    // Get average depth from centroid values
1074    avg_h = ((double *) h -> data)[k];
1075
1076    // Compute bed slope
1077    k6 = 6*k;  // base index
1078 
1079    x0 = ((double*) x -> data)[k6 + 0];
1080    y0 = ((double*) x -> data)[k6 + 1];
1081    x1 = ((double*) x -> data)[k6 + 2];
1082    y1 = ((double*) x -> data)[k6 + 3];
1083    x2 = ((double*) x -> data)[k6 + 4];
1084    y2 = ((double*) x -> data)[k6 + 5];
1085
1086
1087    _gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2, &zx, &zy);
1088
1089    // Update momentum
1090    ((double*) xmom -> data)[k] += -g*zx*avg_h;
1091    ((double*) ymom -> data)[k] += -g*zy*avg_h;
1092  }
1093
1094  return Py_BuildValue("");
1095}
1096
1097
1098PyObject *manning_friction_new(PyObject *self, PyObject *args) {
1099  //
1100  // manning_friction_new(g, eps, x, h, uh, vh, z, eta, xmom_update, ymom_update)
1101  //
1102
1103
1104  PyArrayObject *x, *w, *z, *uh, *vh, *eta, *xmom, *ymom;
1105  int N;
1106  double g, eps;
1107
1108  if (!PyArg_ParseTuple(args, "ddOOOOOOOO",
1109                        &g, &eps, &x, &w, &uh, &vh, &z,  &eta, &xmom, &ymom)) {
1110    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: manning_friction_new could not parse input arguments");
1111    return NULL;
1112  }
1113
1114  // check that numpy array objects arrays are C contiguous memory
1115  CHECK_C_CONTIG(x);
1116  CHECK_C_CONTIG(w);
1117  CHECK_C_CONTIG(z);
1118  CHECK_C_CONTIG(uh);
1119  CHECK_C_CONTIG(vh);
1120  CHECK_C_CONTIG(eta);
1121  CHECK_C_CONTIG(xmom);
1122  CHECK_C_CONTIG(ymom);
1123
1124  N = w -> dimensions[0];
1125
1126  _manning_friction_new(g, eps, N,
1127                        (double*) x -> data,
1128                        (double*) w -> data,
1129                        (double*) z -> data,
1130                        (double*) uh -> data,
1131                        (double*) vh -> data,
1132                        (double*) eta -> data,
1133                        (double*) xmom -> data,
1134                        (double*) ymom -> data);
1135
1136  return Py_BuildValue("");
1137}
1138
1139
1140PyObject *manning_friction_old(PyObject *self, PyObject *args) {
1141  //
1142  // manning_friction_old(g, eps, h, uh, vh, z, eta, xmom_update, ymom_update)
1143  //
1144
1145
1146  PyArrayObject *w, *z, *uh, *vh, *eta, *xmom, *ymom;
1147  int N;
1148  double g, eps;
1149
1150  if (!PyArg_ParseTuple(args, "ddOOOOOOO",
1151                        &g, &eps, &w, &uh, &vh, &z,  &eta, &xmom, &ymom)) {
1152    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: manning_friction_old could not parse input arguments");
1153    return NULL;
1154  }
1155
1156  // check that numpy array objects arrays are C contiguous memory
1157  CHECK_C_CONTIG(w);
1158  CHECK_C_CONTIG(z);
1159  CHECK_C_CONTIG(uh);
1160  CHECK_C_CONTIG(vh);
1161  CHECK_C_CONTIG(eta);
1162  CHECK_C_CONTIG(xmom);
1163  CHECK_C_CONTIG(ymom);
1164
1165  N = w -> dimensions[0];
1166
1167  _manning_friction_old(g, eps, N,
1168                        (double*) w -> data,
1169                        (double*) z -> data,
1170                        (double*) uh -> data,
1171                        (double*) vh -> data,
1172                        (double*) eta -> data,
1173                        (double*) xmom -> data,
1174                        (double*) ymom -> data);
1175
1176  return Py_BuildValue("");
1177}
1178
1179
1180/*
1181PyObject *manning_friction_explicit(PyObject *self, PyObject *args) {
1182  //
1183  // manning_friction_explicit(g, eps, h, uh, vh, eta, xmom_update, ymom_update)
1184  //
1185
1186
1187  PyArrayObject *w, *z, *uh, *vh, *eta, *xmom, *ymom;
1188  int N;
1189  double g, eps;
1190
1191  if (!PyArg_ParseTuple(args, "ddOOOOOOO",
1192            &g, &eps, &w, &z, &uh, &vh, &eta,
1193            &xmom, &ymom))
1194    return NULL;
1195
1196  // check that numpy array objects arrays are C contiguous memory
1197  CHECK_C_CONTIG(w);
1198  CHECK_C_CONTIG(z);
1199  CHECK_C_CONTIG(uh);
1200  CHECK_C_CONTIG(vh);
1201  CHECK_C_CONTIG(eta);
1202  CHECK_C_CONTIG(xmom);
1203  CHECK_C_CONTIG(ymom);
1204
1205  N = w -> dimensions[0];
1206  _manning_friction_explicit(g, eps, N,
1207            (double*) w -> data,
1208            (double*) z -> data,
1209            (double*) uh -> data,
1210            (double*) vh -> data,
1211            (double*) eta -> data,
1212            (double*) xmom -> data,
1213            (double*) ymom -> data);
1214
1215  return Py_BuildValue("");
1216}
1217*/
1218
1219
1220
1221// Computational routine
1222int _extrapolate_second_order_sw(int number_of_elements,
1223                                 double epsilon,
1224                                 double minimum_allowed_height,
1225                                 double beta_w,
1226                                 double beta_w_dry,
1227                                 double beta_uh,
1228                                 double beta_uh_dry,
1229                                 double beta_vh,
1230                                 double beta_vh_dry,
1231                                 long* surrogate_neighbours,
1232                                 long* number_of_boundaries,
1233                                 double* centroid_coordinates,
1234                                 double* stage_centroid_values,
1235                                 double* xmom_centroid_values,
1236                                 double* ymom_centroid_values,
1237                                 double* elevation_centroid_values,
1238                                 double* vertex_coordinates,
1239                                 double* stage_vertex_values,
1240                                 double* xmom_vertex_values,
1241                                 double* ymom_vertex_values,
1242                                 double* elevation_vertex_values,
1243                                 int optimise_dry_cells) {
1244                 
1245                 
1246
1247  // Local variables
1248  double a, b; // Gradient vector used to calculate vertex values from centroids
1249  int k, k0, k1, k2, k3, k6, coord_index, i;
1250  double x, y, x0, y0, x1, y1, x2, y2, xv0, yv0, xv1, yv1, xv2, yv2; // Vertices of the auxiliary triangle
1251  double dx1, dx2, dy1, dy2, dxv0, dxv1, dxv2, dyv0, dyv1, dyv2, dq0, dq1, dq2, area2, inv_area2;
1252  double dqv[3], qmin, qmax, hmin, hmax;
1253  double hc, h0, h1, h2, beta_tmp, hfactor;
1254
1255   
1256  for (k = 0; k < number_of_elements; k++) 
1257  {
1258    k3=k*3;
1259    k6=k*6;
1260
1261    if (number_of_boundaries[k]==3)
1262    {
1263      // No neighbours, set gradient on the triangle to zero
1264     
1265      stage_vertex_values[k3]   = stage_centroid_values[k];
1266      stage_vertex_values[k3+1] = stage_centroid_values[k];
1267      stage_vertex_values[k3+2] = stage_centroid_values[k];
1268      xmom_vertex_values[k3]    = xmom_centroid_values[k];
1269      xmom_vertex_values[k3+1]  = xmom_centroid_values[k];
1270      xmom_vertex_values[k3+2]  = xmom_centroid_values[k];
1271      ymom_vertex_values[k3]    = ymom_centroid_values[k];
1272      ymom_vertex_values[k3+1]  = ymom_centroid_values[k];
1273      ymom_vertex_values[k3+2]  = ymom_centroid_values[k];
1274     
1275      continue;
1276    }
1277    else 
1278    {
1279      // Triangle k has one or more neighbours.
1280      // Get centroid and vertex coordinates of the triangle
1281     
1282      // Get the vertex coordinates
1283      xv0 = vertex_coordinates[k6];   
1284      yv0 = vertex_coordinates[k6+1];
1285      xv1 = vertex_coordinates[k6+2]; 
1286      yv1 = vertex_coordinates[k6+3];
1287      xv2 = vertex_coordinates[k6+4]; 
1288      yv2 = vertex_coordinates[k6+5];
1289     
1290      // Get the centroid coordinates
1291      coord_index = 2*k;
1292      x = centroid_coordinates[coord_index];
1293      y = centroid_coordinates[coord_index+1];
1294     
1295      // Store x- and y- differentials for the vertices of
1296      // triangle k relative to the centroid
1297      dxv0 = xv0 - x; 
1298      dxv1 = xv1 - x; 
1299      dxv2 = xv2 - x;
1300      dyv0 = yv0 - y; 
1301      dyv1 = yv1 - y; 
1302      dyv2 = yv2 - y;
1303    }
1304
1305
1306           
1307    if (number_of_boundaries[k]<=1)
1308    {
1309      //==============================================
1310      // Number of boundaries <= 1
1311      //==============================================   
1312   
1313   
1314      // If no boundaries, auxiliary triangle is formed
1315      // from the centroids of the three neighbours
1316      // If one boundary, auxiliary triangle is formed
1317      // from this centroid and its two neighbours
1318     
1319      k0 = surrogate_neighbours[k3];
1320      k1 = surrogate_neighbours[k3 + 1];
1321      k2 = surrogate_neighbours[k3 + 2];
1322     
1323      // Get the auxiliary triangle's vertex coordinates
1324      // (really the centroids of neighbouring triangles)
1325      coord_index = 2*k0;
1326      x0 = centroid_coordinates[coord_index];
1327      y0 = centroid_coordinates[coord_index+1];
1328     
1329      coord_index = 2*k1;
1330      x1 = centroid_coordinates[coord_index];
1331      y1 = centroid_coordinates[coord_index+1];
1332     
1333      coord_index = 2*k2;
1334      x2 = centroid_coordinates[coord_index];
1335      y2 = centroid_coordinates[coord_index+1];
1336     
1337      // Store x- and y- differentials for the vertices
1338      // of the auxiliary triangle
1339      dx1 = x1 - x0; 
1340      dx2 = x2 - x0;
1341      dy1 = y1 - y0; 
1342      dy2 = y2 - y0;
1343     
1344      // Calculate 2*area of the auxiliary triangle
1345      // The triangle is guaranteed to be counter-clockwise     
1346      area2 = dy2*dx1 - dy1*dx2; 
1347     
1348      // If the mesh is 'weird' near the boundary,
1349      // the triangle might be flat or clockwise:
1350      if (area2 <= 0) 
1351      {
1352        PyErr_SetString(PyExc_RuntimeError, 
1353                        "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.");
1354   
1355        return -1;
1356      } 
1357     
1358      // Calculate heights of neighbouring cells
1359      hc = stage_centroid_values[k]  - elevation_centroid_values[k];
1360      h0 = stage_centroid_values[k0] - elevation_centroid_values[k0];
1361      h1 = stage_centroid_values[k1] - elevation_centroid_values[k1];
1362      h2 = stage_centroid_values[k2] - elevation_centroid_values[k2];
1363      hmin = min(min(h0, min(h1, h2)), hc);
1364      //hfactor = hc/(hc + 1.0);
1365
1366      hfactor = 0.0;
1367      if (hmin > 0.001) 
1368      {
1369        hfactor = (hmin - 0.001)/(hmin + 0.004);
1370      }
1371     
1372      if (optimise_dry_cells) 
1373      {     
1374        // Check if linear reconstruction is necessary for triangle k
1375        // This check will exclude dry cells.
1376
1377        hmax = max(h0, max(h1, h2));     
1378        if (hmax < epsilon) 
1379        {
1380            continue;
1381        }
1382      }
1383   
1384      //-----------------------------------
1385      // stage
1386      //-----------------------------------     
1387     
1388      // Calculate the difference between vertex 0 of the auxiliary
1389      // triangle and the centroid of triangle k
1390      dq0 = stage_centroid_values[k0] - stage_centroid_values[k];
1391     
1392      // Calculate differentials between the vertices
1393      // of the auxiliary triangle (centroids of neighbouring triangles)
1394      dq1 = stage_centroid_values[k1] - stage_centroid_values[k0];
1395      dq2 = stage_centroid_values[k2] - stage_centroid_values[k0];
1396     
1397          inv_area2 = 1.0/area2;
1398      // Calculate the gradient of stage on the auxiliary triangle
1399      a = dy2*dq1 - dy1*dq2;
1400      a *= inv_area2;
1401      b = dx1*dq2 - dx2*dq1;
1402      b *= inv_area2;
1403     
1404      // Calculate provisional jumps in stage from the centroid
1405      // of triangle k to its vertices, to be limited
1406      dqv[0] = a*dxv0 + b*dyv0;
1407      dqv[1] = a*dxv1 + b*dyv1;
1408      dqv[2] = a*dxv2 + b*dyv2;
1409     
1410      // Now we want to find min and max of the centroid and the
1411      // vertices of the auxiliary triangle and compute jumps
1412      // from the centroid to the min and max
1413      find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
1414     
1415      // Playing with dry wet interface
1416      //hmin = qmin;
1417      //beta_tmp = beta_w_dry;
1418      //if (hmin>minimum_allowed_height)
1419      beta_tmp = beta_w_dry + (beta_w - beta_w_dry) * hfactor;
1420   
1421      //printf("min_alled_height = %f\n",minimum_allowed_height);
1422      //printf("hmin = %f\n",hmin);
1423      //printf("beta_w = %f\n",beta_w);
1424      //printf("beta_tmp = %f\n",beta_tmp);
1425      // Limit the gradient
1426      limit_gradient(dqv, qmin, qmax, beta_tmp); 
1427     
1428      //for (i=0;i<3;i++)
1429      stage_vertex_values[k3+0] = stage_centroid_values[k] + dqv[0];
1430          stage_vertex_values[k3+1] = stage_centroid_values[k] + dqv[1];
1431          stage_vertex_values[k3+2] = stage_centroid_values[k] + dqv[2];
1432     
1433     
1434      //-----------------------------------
1435      // xmomentum
1436      //-----------------------------------           
1437
1438      // Calculate the difference between vertex 0 of the auxiliary
1439      // triangle and the centroid of triangle k     
1440      dq0 = xmom_centroid_values[k0] - xmom_centroid_values[k];
1441     
1442      // Calculate differentials between the vertices
1443      // of the auxiliary triangle
1444      dq1 = xmom_centroid_values[k1] - xmom_centroid_values[k0];
1445      dq2 = xmom_centroid_values[k2] - xmom_centroid_values[k0];
1446     
1447      // Calculate the gradient of xmom on the auxiliary triangle
1448      a = dy2*dq1 - dy1*dq2;
1449      a *= inv_area2;
1450      b = dx1*dq2 - dx2*dq1;
1451      b *= inv_area2;
1452     
1453      // Calculate provisional jumps in stage from the centroid
1454      // of triangle k to its vertices, to be limited     
1455      dqv[0] = a*dxv0+b*dyv0;
1456      dqv[1] = a*dxv1+b*dyv1;
1457      dqv[2] = a*dxv2+b*dyv2;
1458     
1459      // Now we want to find min and max of the centroid and the
1460      // vertices of the auxiliary triangle and compute jumps
1461      // from the centroid to the min and max
1462      find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
1463      //beta_tmp = beta_uh;
1464      //if (hmin<minimum_allowed_height)
1465      //beta_tmp = beta_uh_dry;
1466      beta_tmp = beta_uh_dry + (beta_uh - beta_uh_dry) * hfactor;
1467
1468      // Limit the gradient
1469      limit_gradient(dqv, qmin, qmax, beta_tmp);
1470
1471      for (i=0; i < 3; i++)
1472      {
1473        xmom_vertex_values[k3+i] = xmom_centroid_values[k] + dqv[i];
1474      }
1475     
1476      //-----------------------------------
1477      // ymomentum
1478      //-----------------------------------                 
1479
1480      // Calculate the difference between vertex 0 of the auxiliary
1481      // triangle and the centroid of triangle k
1482      dq0 = ymom_centroid_values[k0] - ymom_centroid_values[k];
1483     
1484      // Calculate differentials between the vertices
1485      // of the auxiliary triangle
1486      dq1 = ymom_centroid_values[k1] - ymom_centroid_values[k0];
1487      dq2 = ymom_centroid_values[k2] - ymom_centroid_values[k0];
1488     
1489      // Calculate the gradient of xmom on the auxiliary triangle
1490      a = dy2*dq1 - dy1*dq2;
1491      a *= inv_area2;
1492      b = dx1*dq2 - dx2*dq1;
1493      b *= inv_area2;
1494     
1495      // Calculate provisional jumps in stage from the centroid
1496      // of triangle k to its vertices, to be limited
1497      dqv[0] = a*dxv0 + b*dyv0;
1498      dqv[1] = a*dxv1 + b*dyv1;
1499      dqv[2] = a*dxv2 + b*dyv2;
1500     
1501      // Now we want to find min and max of the centroid and the
1502      // vertices of the auxiliary triangle and compute jumps
1503      // from the centroid to the min and max
1504      find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
1505     
1506      //beta_tmp = beta_vh;
1507      //
1508      //if (hmin<minimum_allowed_height)
1509      //beta_tmp = beta_vh_dry;
1510      beta_tmp = beta_vh_dry + (beta_vh - beta_vh_dry) * hfactor;   
1511
1512      // Limit the gradient
1513      limit_gradient(dqv, qmin, qmax, beta_tmp);
1514     
1515      for (i=0;i<3;i++)
1516      {
1517        ymom_vertex_values[k3 + i] = ymom_centroid_values[k] + dqv[i];
1518      }
1519    } // End number_of_boundaries <=1
1520    else
1521    {
1522
1523      //==============================================
1524      // Number of boundaries == 2
1525      //==============================================       
1526       
1527      // One internal neighbour and gradient is in direction of the neighbour's centroid
1528     
1529      // Find the only internal neighbour (k1?)
1530      for (k2 = k3; k2 < k3 + 3; k2++)
1531      {
1532      // Find internal neighbour of triangle k     
1533      // k2 indexes the edges of triangle k   
1534     
1535          if (surrogate_neighbours[k2] != k)
1536          {
1537             break;
1538          }
1539      }
1540     
1541      if ((k2 == k3 + 3)) 
1542      {
1543        // If we didn't find an internal neighbour
1544        PyErr_SetString(PyExc_RuntimeError, 
1545                        "shallow_water_ext.c: Internal neighbour not found");     
1546        return -1;
1547      }
1548     
1549      k1 = surrogate_neighbours[k2];
1550     
1551      // The coordinates of the triangle are already (x,y).
1552      // Get centroid of the neighbour (x1,y1)
1553      coord_index = 2*k1;
1554      x1 = centroid_coordinates[coord_index];
1555      y1 = centroid_coordinates[coord_index + 1];
1556     
1557      // Compute x- and y- distances between the centroid of
1558      // triangle k and that of its neighbour
1559      dx1 = x1 - x; 
1560      dy1 = y1 - y;
1561     
1562      // Set area2 as the square of the distance
1563      area2 = dx1*dx1 + dy1*dy1;
1564     
1565      // Set dx2=(x1-x0)/((x1-x0)^2+(y1-y0)^2)
1566      // and dy2=(y1-y0)/((x1-x0)^2+(y1-y0)^2) which
1567      // respectively correspond to the x- and y- gradients
1568      // of the conserved quantities
1569      dx2 = 1.0/area2;
1570      dy2 = dx2*dy1;
1571      dx2 *= dx1;
1572     
1573     
1574      //-----------------------------------
1575      // stage
1576      //-----------------------------------           
1577
1578      // Compute differentials
1579      dq1 = stage_centroid_values[k1] - stage_centroid_values[k];
1580     
1581      // Calculate the gradient between the centroid of triangle k
1582      // and that of its neighbour
1583      a = dq1*dx2;
1584      b = dq1*dy2;
1585     
1586      // Calculate provisional vertex jumps, to be limited
1587      dqv[0] = a*dxv0 + b*dyv0;
1588      dqv[1] = a*dxv1 + b*dyv1;
1589      dqv[2] = a*dxv2 + b*dyv2;
1590     
1591      // Now limit the jumps
1592      if (dq1>=0.0)
1593      {
1594        qmin=0.0;
1595        qmax=dq1;
1596      }
1597      else
1598      {
1599        qmin = dq1;
1600        qmax = 0.0;
1601      }
1602     
1603      // Limit the gradient
1604      limit_gradient(dqv, qmin, qmax, beta_w);
1605     
1606      //for (i=0; i < 3; i++)
1607      //{
1608      stage_vertex_values[k3] = stage_centroid_values[k] + dqv[0];
1609      stage_vertex_values[k3 + 1] = stage_centroid_values[k] + dqv[1];
1610      stage_vertex_values[k3 + 2] = stage_centroid_values[k] + dqv[2];
1611      //}
1612     
1613      //-----------------------------------
1614      // xmomentum
1615      //-----------------------------------                       
1616     
1617      // Compute differentials
1618      dq1 = xmom_centroid_values[k1] - xmom_centroid_values[k];
1619     
1620      // Calculate the gradient between the centroid of triangle k
1621      // and that of its neighbour
1622      a = dq1*dx2;
1623      b = dq1*dy2;
1624     
1625      // Calculate provisional vertex jumps, to be limited
1626      dqv[0] = a*dxv0+b*dyv0;
1627      dqv[1] = a*dxv1+b*dyv1;
1628      dqv[2] = a*dxv2+b*dyv2;
1629     
1630      // Now limit the jumps
1631      if (dq1 >= 0.0)
1632      {
1633        qmin = 0.0;
1634        qmax = dq1;
1635      }
1636      else
1637      {
1638        qmin = dq1;
1639        qmax = 0.0;
1640      }
1641     
1642      // Limit the gradient     
1643      limit_gradient(dqv, qmin, qmax, beta_w);
1644     
1645      //for (i=0;i<3;i++)
1646      //xmom_vertex_values[k3] = xmom_centroid_values[k] + dqv[0];
1647      //xmom_vertex_values[k3 + 1] = xmom_centroid_values[k] + dqv[1];
1648      //xmom_vertex_values[k3 + 2] = xmom_centroid_values[k] + dqv[2];
1649     
1650      for (i = 0; i < 3;i++)
1651      {
1652          xmom_vertex_values[k3 + i] = xmom_centroid_values[k] + dqv[i];
1653      }
1654     
1655      //-----------------------------------
1656      // ymomentum
1657      //-----------------------------------                       
1658
1659      // Compute differentials
1660      dq1 = ymom_centroid_values[k1] - ymom_centroid_values[k];
1661     
1662      // Calculate the gradient between the centroid of triangle k
1663      // and that of its neighbour
1664      a = dq1*dx2;
1665      b = dq1*dy2;
1666     
1667      // Calculate provisional vertex jumps, to be limited
1668      dqv[0] = a*dxv0 + b*dyv0;
1669      dqv[1] = a*dxv1 + b*dyv1;
1670      dqv[2] = a*dxv2 + b*dyv2;
1671     
1672      // Now limit the jumps
1673      if (dq1>=0.0)
1674      {
1675        qmin = 0.0;
1676        qmax = dq1;
1677      } 
1678      else
1679      {
1680        qmin = dq1;
1681        qmax = 0.0;
1682      }
1683     
1684      // Limit the gradient
1685      limit_gradient(dqv, qmin, qmax, beta_w);
1686     
1687      //for (i=0;i<3;i++)
1688      //ymom_vertex_values[k3] = ymom_centroid_values[k] + dqv[0];
1689      //ymom_vertex_values[k3 + 1] = ymom_centroid_values[k] + dqv[1];
1690      //ymom_vertex_values[k3 + 2] = ymom_centroid_values[k] + dqv[2];
1691
1692for (i=0;i<3;i++)
1693        {
1694ymom_vertex_values[k3 + i] = ymom_centroid_values[k] + dqv[i];
1695        }
1696//ymom_vertex_values[k3] = ymom_centroid_values[k] + dqv[0];
1697//ymom_vertex_values[k3 + 1] = ymom_centroid_values[k] + dqv[1];
1698//ymom_vertex_values[k3 + 2] = ymom_centroid_values[k] + dqv[2];
1699    } // else [number_of_boundaries==2]
1700  } // for k=0 to number_of_elements-1
1701 
1702  return 0;
1703}           
1704
1705
1706PyObject *extrapolate_second_order_sw(PyObject *self, PyObject *args) {
1707  /*Compute the vertex values based on a linear reconstruction
1708    on each triangle
1709   
1710    These values are calculated as follows:
1711    1) For each triangle not adjacent to a boundary, we consider the
1712       auxiliary triangle formed by the centroids of its three
1713       neighbours.
1714    2) For each conserved quantity, we integrate around the auxiliary
1715       triangle's boundary the product of the quantity and the outward
1716       normal vector. Dividing by the triangle area gives (a,b), the
1717       average of the vector (q_x,q_y) on the auxiliary triangle.
1718       We suppose that the linear reconstruction on the original
1719       triangle has gradient (a,b).
1720    3) Provisional vertex jumps dqv[0,1,2] are computed and these are
1721       then limited by calling the functions find_qmin_and_qmax and
1722       limit_gradient
1723
1724    Python call:
1725    extrapolate_second_order_sw(domain.surrogate_neighbours,
1726                                domain.number_of_boundaries
1727                                domain.centroid_coordinates,
1728                                Stage.centroid_values
1729                                Xmom.centroid_values
1730                                Ymom.centroid_values
1731                                domain.vertex_coordinates,
1732                                Stage.vertex_values,
1733                                Xmom.vertex_values,
1734                                Ymom.vertex_values)
1735
1736    Post conditions:
1737            The vertices of each triangle have values from a
1738        limited linear reconstruction
1739        based on centroid values
1740
1741  */
1742  PyArrayObject *surrogate_neighbours,
1743    *number_of_boundaries,
1744    *centroid_coordinates,
1745    *stage_centroid_values,
1746    *xmom_centroid_values,
1747    *ymom_centroid_values,
1748    *elevation_centroid_values,
1749    *vertex_coordinates,
1750    *stage_vertex_values,
1751    *xmom_vertex_values,
1752    *ymom_vertex_values,
1753    *elevation_vertex_values;
1754 
1755  PyObject *domain;
1756
1757 
1758  double beta_w, beta_w_dry, beta_uh, beta_uh_dry, beta_vh, beta_vh_dry;   
1759  double minimum_allowed_height, epsilon;
1760  int optimise_dry_cells, number_of_elements, e;
1761 
1762  // Provisional jumps from centroids to v'tices and safety factor re limiting
1763  // by which these jumps are limited
1764  // Convert Python arguments to C
1765  if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOi",
1766            &domain,
1767            &surrogate_neighbours,
1768            &number_of_boundaries,
1769            &centroid_coordinates,
1770            &stage_centroid_values,
1771            &xmom_centroid_values,
1772            &ymom_centroid_values,
1773            &elevation_centroid_values,
1774            &vertex_coordinates,
1775            &stage_vertex_values,
1776            &xmom_vertex_values,
1777            &ymom_vertex_values,
1778            &elevation_vertex_values,
1779            &optimise_dry_cells)) {         
1780           
1781    PyErr_SetString(PyExc_RuntimeError, 
1782            "Input arguments to extrapolate_second_order_sw failed");
1783    return NULL;
1784  }
1785
1786  // check that numpy array objects arrays are C contiguous memory
1787  CHECK_C_CONTIG(surrogate_neighbours);
1788  CHECK_C_CONTIG(number_of_boundaries);
1789  CHECK_C_CONTIG(centroid_coordinates);
1790  CHECK_C_CONTIG(stage_centroid_values);
1791  CHECK_C_CONTIG(xmom_centroid_values);
1792  CHECK_C_CONTIG(ymom_centroid_values);
1793  CHECK_C_CONTIG(elevation_centroid_values);
1794  CHECK_C_CONTIG(vertex_coordinates);
1795  CHECK_C_CONTIG(stage_vertex_values);
1796  CHECK_C_CONTIG(xmom_vertex_values);
1797  CHECK_C_CONTIG(ymom_vertex_values);
1798  CHECK_C_CONTIG(elevation_vertex_values);
1799 
1800  // Get the safety factor beta_w, set in the config.py file.
1801  // This is used in the limiting process
1802 
1803
1804  beta_w                 = get_python_double(domain,"beta_w");
1805  beta_w_dry             = get_python_double(domain,"beta_w_dry");
1806  beta_uh                = get_python_double(domain,"beta_uh");
1807  beta_uh_dry            = get_python_double(domain,"beta_uh_dry");
1808  beta_vh                = get_python_double(domain,"beta_vh");
1809  beta_vh_dry            = get_python_double(domain,"beta_vh_dry"); 
1810
1811  minimum_allowed_height = get_python_double(domain,"minimum_allowed_height");
1812  epsilon                = get_python_double(domain,"epsilon");
1813
1814  number_of_elements = stage_centroid_values -> dimensions[0]; 
1815
1816  // Call underlying computational routine
1817  e = _extrapolate_second_order_sw(number_of_elements,
1818                   epsilon,
1819                   minimum_allowed_height,
1820                   beta_w,
1821                   beta_w_dry,
1822                   beta_uh,
1823                   beta_uh_dry,
1824                   beta_vh,
1825                   beta_vh_dry,
1826                   (long*) surrogate_neighbours -> data,
1827                   (long*) number_of_boundaries -> data,
1828                   (double*) centroid_coordinates -> data,
1829                   (double*) stage_centroid_values -> data,
1830                   (double*) xmom_centroid_values -> data,
1831                   (double*) ymom_centroid_values -> data,
1832                   (double*) elevation_centroid_values -> data,
1833                   (double*) vertex_coordinates -> data,
1834                   (double*) stage_vertex_values -> data,
1835                   (double*) xmom_vertex_values -> data,
1836                   (double*) ymom_vertex_values -> data,
1837                   (double*) elevation_vertex_values -> data,
1838                   optimise_dry_cells);
1839  if (e == -1) {
1840    // Use error string set inside computational routine
1841    return NULL;
1842  }                   
1843 
1844 
1845  return Py_BuildValue("");
1846 
1847}// extrapolate_second-order_sw
1848
1849
1850
1851
1852PyObject *rotate(PyObject *self, PyObject *args, PyObject *kwargs) {
1853  //
1854  // r = rotate(q, normal, direction=1)
1855  //
1856  // Where q is assumed to be a Float numeric array of length 3 and
1857  // normal a Float numeric array of length 2.
1858
1859  // FIXME(Ole): I don't think this is used anymore
1860
1861  PyObject *Q, *Normal;
1862  PyArrayObject *q, *r, *normal;
1863
1864  static char *argnames[] = {"q", "normal", "direction", NULL};
1865  int dimensions[1], i, direction=1;
1866  double n1, n2;
1867
1868  // Convert Python arguments to C
1869  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OO|i", argnames,
1870                   &Q, &Normal, &direction)) {
1871    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: rotate could not parse input arguments");
1872    return NULL;
1873  } 
1874
1875  // Input checks (convert sequences into numeric arrays)
1876  q = (PyArrayObject *)
1877    PyArray_ContiguousFromObject(Q, PyArray_DOUBLE, 0, 0);
1878  normal = (PyArrayObject *)
1879    PyArray_ContiguousFromObject(Normal, PyArray_DOUBLE, 0, 0);
1880
1881
1882  if (normal -> dimensions[0] != 2) {
1883    PyErr_SetString(PyExc_RuntimeError, "Normal vector must have 2 components");
1884    return NULL;
1885  }
1886
1887  // Allocate space for return vector r (don't DECREF)
1888  dimensions[0] = 3;
1889  r = (PyArrayObject *) anuga_FromDims(1, dimensions, PyArray_DOUBLE);
1890
1891  // Copy
1892  for (i=0; i<3; i++) {
1893    ((double *) (r -> data))[i] = ((double *) (q -> data))[i];
1894  }
1895
1896  // Get normal and direction
1897  n1 = ((double *) normal -> data)[0];
1898  n2 = ((double *) normal -> data)[1];
1899  if (direction == -1) n2 = -n2;
1900
1901  // Rotate
1902  _rotate((double *) r -> data, n1, n2);
1903
1904  // Release numeric arrays
1905  Py_DECREF(q);
1906  Py_DECREF(normal);
1907
1908  // Return result using PyArray to avoid memory leak
1909  return PyArray_Return(r);
1910}
1911
1912
1913// Computational function for flux computation
1914double _compute_fluxes_central(int number_of_elements,
1915                               double timestep,
1916                               double epsilon,
1917                               double H0,
1918                               double g,
1919                               long* neighbours,
1920                               long* neighbour_edges,
1921                               double* normals,
1922                               double* edgelengths, 
1923                               double* radii, 
1924                               double* areas,
1925                               long* tri_full_flag,
1926                               double* stage_edge_values,
1927                               double* xmom_edge_values,
1928                               double* ymom_edge_values,
1929                               double* bed_edge_values,
1930                               double* stage_boundary_values,
1931                               double* xmom_boundary_values,
1932                               double* ymom_boundary_values,
1933                               double* stage_explicit_update,
1934                               double* xmom_explicit_update,
1935                               double* ymom_explicit_update,
1936                               long* already_computed_flux,
1937                               double* max_speed_array,
1938                               int optimise_dry_cells) 
1939{                   
1940  // Local variables
1941  double max_speed, length, inv_area, zl, zr;
1942  double h0 = H0*H0; // This ensures a good balance when h approaches H0.
1943
1944  double limiting_threshold = 10*H0; // Avoid applying limiter below this
1945                                     // threshold for performance reasons.
1946                                     // See ANUGA manual under flux limiting
1947  int k, i, m, n;
1948  int ki, nm=0, ki2; // Index shorthands
1949 
1950 
1951  // Workspace (making them static actually made function slightly slower (Ole)) 
1952  double ql[3], qr[3], edgeflux[3]; // Work array for summing up fluxes
1953
1954  static long call = 1; // Static local variable flagging already computed flux
1955 
1956  // Start computation
1957  call++; // Flag 'id' of flux calculation for this timestep
1958
1959  // Set explicit_update to zero for all conserved_quantities.
1960  // This assumes compute_fluxes called before forcing terms
1961  memset((char*) stage_explicit_update, 0, number_of_elements*sizeof(double));
1962  memset((char*) xmom_explicit_update, 0, number_of_elements*sizeof(double));
1963  memset((char*) ymom_explicit_update, 0, number_of_elements*sizeof(double));   
1964 
1965  // For all triangles
1966  for (k = 0; k < number_of_elements; k++) 
1967  { 
1968    // Loop through neighbours and compute edge flux for each 
1969    for (i = 0; i < 3; i++) 
1970    {
1971      ki = k*3 + i; // Linear index to edge i of triangle k
1972     
1973      if (already_computed_flux[ki] == call)
1974      {
1975        // We've already computed the flux across this edge
1976        continue;
1977      }
1978   
1979      // Get left hand side values from triangle k, edge i
1980      ql[0] = stage_edge_values[ki];
1981      ql[1] = xmom_edge_values[ki];
1982      ql[2] = ymom_edge_values[ki];
1983      zl = bed_edge_values[ki];
1984
1985      // Get right hand side values either from neighbouring triangle
1986      // or from boundary array (Quantities at neighbour on nearest face).
1987      n = neighbours[ki];
1988      if (n < 0) 
1989      {
1990        // Neighbour is a boundary condition
1991        m = -n - 1; // Convert negative flag to boundary index
1992   
1993        qr[0] = stage_boundary_values[m];
1994        qr[1] = xmom_boundary_values[m];
1995        qr[2] = ymom_boundary_values[m];
1996        zr = zl; // Extend bed elevation to boundary
1997      } 
1998      else 
1999      {
2000        // Neighbour is a real triangle
2001        m = neighbour_edges[ki];
2002        nm = n*3 + m; // Linear index (triangle n, edge m)
2003       
2004        qr[0] = stage_edge_values[nm];
2005        qr[1] = xmom_edge_values[nm];
2006        qr[2] = ymom_edge_values[nm];
2007        zr = bed_edge_values[nm];
2008      }
2009     
2010      // Now we have values for this edge - both from left and right side.
2011     
2012      if (optimise_dry_cells) 
2013      {     
2014        // Check if flux calculation is necessary across this edge
2015        // This check will exclude dry cells.
2016        // This will also optimise cases where zl != zr as
2017        // long as both are dry
2018
2019        if (fabs(ql[0] - zl) < epsilon && 
2020            fabs(qr[0] - zr) < epsilon) 
2021        {
2022          // Cell boundary is dry
2023         
2024          already_computed_flux[ki] = call; // #k Done 
2025          if (n >= 0)
2026          {
2027            already_computed_flux[nm] = call; // #n Done
2028          }
2029       
2030          max_speed = 0.0;
2031          continue;
2032        }
2033      }
2034     
2035      // Outward pointing normal vector (domain.normals[k, 2*i:2*i+2])
2036      ki2 = 2*ki; //k*6 + i*2
2037
2038      // Edge flux computation (triangle k, edge i)
2039      _flux_function_central(ql, qr, zl, zr,
2040                             normals[ki2], normals[ki2+1],
2041                             epsilon, h0, limiting_threshold, g,
2042                             edgeflux, &max_speed);
2043     
2044     
2045      // Multiply edgeflux by edgelength
2046      length = edgelengths[ki];
2047      edgeflux[0] *= length;           
2048      edgeflux[1] *= length;           
2049      edgeflux[2] *= length;                       
2050     
2051     
2052      // Update triangle k with flux from edge i
2053      stage_explicit_update[k] -= edgeflux[0];
2054      xmom_explicit_update[k] -= edgeflux[1];
2055      ymom_explicit_update[k] -= edgeflux[2];
2056     
2057      already_computed_flux[ki] = call; // #k Done
2058     
2059     
2060      // Update neighbour n with same flux but reversed sign
2061      if (n >= 0) 
2062      {
2063        stage_explicit_update[n] += edgeflux[0];
2064        xmom_explicit_update[n] += edgeflux[1];
2065        ymom_explicit_update[n] += edgeflux[2];
2066       
2067        already_computed_flux[nm] = call; // #n Done
2068      }
2069
2070      // Update timestep based on edge i and possibly neighbour n
2071      if (tri_full_flag[k] == 1) 
2072      {
2073        if (max_speed > epsilon) 
2074        {
2075          // Apply CFL condition for triangles joining this edge (triangle k and triangle n)
2076         
2077          // CFL for triangle k
2078          timestep = min(timestep, radii[k]/max_speed);
2079         
2080          if (n >= 0)
2081          {
2082            // Apply CFL condition for neigbour n (which is on the ith edge of triangle k)
2083            timestep = min(timestep, radii[n]/max_speed);
2084          }
2085
2086          // Ted Rigby's suggested less conservative version
2087          //if (n>=0) {         
2088          //  timestep = min(timestep, (radii[k]+radii[n])/max_speed);
2089          //} else {
2090          //  timestep = min(timestep, radii[k]/max_speed);
2091          // }
2092        }
2093      }
2094     
2095    } // End edge i (and neighbour n)
2096   
2097   
2098    // Normalise triangle k by area and store for when all conserved
2099    // quantities get updated
2100    inv_area = 1.0/areas[k];
2101    stage_explicit_update[k] *= inv_area;
2102    xmom_explicit_update[k] *= inv_area;
2103    ymom_explicit_update[k] *= inv_area;
2104   
2105   
2106    // Keep track of maximal speeds
2107    max_speed_array[k] = max_speed;
2108   
2109  } // End triangle k
2110             
2111  return timestep;
2112}
2113
2114//=========================================================================
2115// Python Glue
2116//=========================================================================
2117
2118PyObject *compute_fluxes_ext_central_new(PyObject *self, PyObject *args) {
2119  /*Compute all fluxes and the timestep suitable for all volumes
2120    in domain.
2121
2122    Compute total flux for each conserved quantity using "flux_function_central"
2123
2124    Fluxes across each edge are scaled by edgelengths and summed up
2125    Resulting flux is then scaled by area and stored in
2126    explicit_update for each of the three conserved quantities
2127    stage, xmomentum and ymomentum
2128
2129    The maximal allowable speed computed by the flux_function for each volume
2130    is converted to a timestep that must not be exceeded. The minimum of
2131    those is computed as the next overall timestep.
2132
2133    Python call:
2134    timestep = compute_fluxes(timestep, domain, stage, xmom, ymom, bed)
2135
2136
2137    Post conditions:
2138      domain.explicit_update is reset to computed flux values
2139      returns timestep which is the largest step satisfying all volumes.
2140
2141
2142  */
2143
2144    PyObject
2145    *domain,
2146    *stage, 
2147    *xmom, 
2148    *ymom, 
2149    *bed;
2150
2151    PyArrayObject
2152    *neighbours, 
2153    *neighbour_edges,
2154    *normals, 
2155    *edgelengths, 
2156    *radii, 
2157    *areas,
2158    *tri_full_flag,
2159    *stage_edge_values,
2160    *xmom_edge_values,
2161    *ymom_edge_values,
2162    *bed_edge_values,
2163    *stage_boundary_values,
2164    *xmom_boundary_values,
2165    *ymom_boundary_values,
2166    *stage_explicit_update,
2167    *xmom_explicit_update,
2168    *ymom_explicit_update,
2169    *already_computed_flux, //Tracks whether the flux across an edge has already been computed
2170    *max_speed_array; //Keeps track of max speeds for each triangle
2171
2172   
2173    double timestep, epsilon, H0, g;
2174    int optimise_dry_cells;
2175   
2176    // Convert Python arguments to C
2177    if (!PyArg_ParseTuple(args, "dOOOO", &timestep, &domain, &stage, &xmom, &ymom, &bed )) {
2178    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
2179    return NULL;
2180    }
2181
2182    epsilon           = get_python_double(domain,"epsilon");
2183    H0                = get_python_double(domain,"H0");
2184    g                 = get_python_double(domain,"g");
2185    optimise_dry_cells = get_python_integer(domain,"optimse_dry_cells");
2186   
2187    neighbours        = get_consecutive_array(domain, "neighbours");
2188    neighbour_edges   = get_consecutive_array(domain, "neighbour_edges"); 
2189    normals           = get_consecutive_array(domain, "normals");
2190    edgelengths       = get_consecutive_array(domain, "edge_lengths");   
2191    radii             = get_consecutive_array(domain, "radii");   
2192    areas             = get_consecutive_array(domain, "areas");   
2193    tri_full_flag     = get_consecutive_array(domain, "tri_full_flag");
2194    already_computed_flux  = get_consecutive_array(domain, "already_computed_flux");
2195    max_speed_array   = get_consecutive_array(domain, "max_speed");
2196   
2197    stage_edge_values = get_consecutive_array(stage, "edge_values");   
2198    xmom_edge_values  = get_consecutive_array(xmom, "edge_values");   
2199    ymom_edge_values  = get_consecutive_array(ymom, "edge_values"); 
2200    bed_edge_values   = get_consecutive_array(bed, "edge_values");   
2201
2202    stage_boundary_values = get_consecutive_array(stage, "boundary_values");   
2203    xmom_boundary_values  = get_consecutive_array(xmom, "boundary_values");   
2204    ymom_boundary_values  = get_consecutive_array(ymom, "boundary_values"); 
2205
2206    stage_explicit_update = get_consecutive_array(stage, "explicit_update");   
2207    xmom_explicit_update  = get_consecutive_array(xmom, "explicit_update");   
2208    ymom_explicit_update  = get_consecutive_array(ymom, "explicit_update"); 
2209
2210
2211  int number_of_elements = stage_edge_values -> dimensions[0];
2212
2213  // Call underlying flux computation routine and update
2214  // the explicit update arrays
2215  timestep = _compute_fluxes_central(number_of_elements,
2216                     timestep,
2217                     epsilon,
2218                     H0,
2219                     g,
2220                     (long*) neighbours -> data,
2221                     (long*) neighbour_edges -> data,
2222                     (double*) normals -> data,
2223                     (double*) edgelengths -> data, 
2224                     (double*) radii -> data, 
2225                     (double*) areas -> data,
2226                     (long*) tri_full_flag -> data,
2227                     (double*) stage_edge_values -> data,
2228                     (double*) xmom_edge_values -> data,
2229                     (double*) ymom_edge_values -> data,
2230                     (double*) bed_edge_values -> data,
2231                     (double*) stage_boundary_values -> data,
2232                     (double*) xmom_boundary_values -> data,
2233                     (double*) ymom_boundary_values -> data,
2234                     (double*) stage_explicit_update -> data,
2235                     (double*) xmom_explicit_update -> data,
2236                     (double*) ymom_explicit_update -> data,
2237                     (long*) already_computed_flux -> data,
2238                     (double*) max_speed_array -> data,
2239                     optimise_dry_cells);
2240
2241  Py_DECREF(neighbours);
2242  Py_DECREF(neighbour_edges);
2243  Py_DECREF(normals);
2244  Py_DECREF(edgelengths);
2245  Py_DECREF(radii);
2246  Py_DECREF(areas);
2247  Py_DECREF(tri_full_flag);
2248  Py_DECREF(already_computed_flux);
2249  Py_DECREF(max_speed_array);
2250  Py_DECREF(stage_edge_values);
2251  Py_DECREF(xmom_edge_values);
2252  Py_DECREF(ymom_edge_values);
2253  Py_DECREF(bed_edge_values);
2254  Py_DECREF(stage_boundary_values);
2255  Py_DECREF(xmom_boundary_values);
2256  Py_DECREF(ymom_boundary_values);
2257  Py_DECREF(stage_explicit_update);
2258  Py_DECREF(xmom_explicit_update);
2259  Py_DECREF(ymom_explicit_update);
2260
2261 
2262  // Return updated flux timestep
2263  return Py_BuildValue("d", timestep);
2264}
2265
2266
2267
2268
2269
2270
2271PyObject *compute_fluxes_ext_central(PyObject *self, PyObject *args) {
2272  /*Compute all fluxes and the timestep suitable for all volumes
2273    in domain.
2274
2275    Compute total flux for each conserved quantity using "flux_function_central"
2276
2277    Fluxes across each edge are scaled by edgelengths and summed up
2278    Resulting flux is then scaled by area and stored in
2279    explicit_update for each of the three conserved quantities
2280    stage, xmomentum and ymomentum
2281
2282    The maximal allowable speed computed by the flux_function for each volume
2283    is converted to a timestep that must not be exceeded. The minimum of
2284    those is computed as the next overall timestep.
2285
2286    Python call:
2287    domain.timestep = compute_fluxes(timestep,
2288                                     domain.epsilon,
2289                     domain.H0,
2290                                     domain.g,
2291                                     domain.neighbours,
2292                                     domain.neighbour_edges,
2293                                     domain.normals,
2294                                     domain.edgelengths,
2295                                     domain.radii,
2296                                     domain.areas,
2297                                     tri_full_flag,
2298                                     Stage.edge_values,
2299                                     Xmom.edge_values,
2300                                     Ymom.edge_values,
2301                                     Bed.edge_values,
2302                                     Stage.boundary_values,
2303                                     Xmom.boundary_values,
2304                                     Ymom.boundary_values,
2305                                     Stage.explicit_update,
2306                                     Xmom.explicit_update,
2307                                     Ymom.explicit_update,
2308                                     already_computed_flux,
2309                     optimise_dry_cells)                     
2310
2311
2312    Post conditions:
2313      domain.explicit_update is reset to computed flux values
2314      domain.timestep is set to the largest step satisfying all volumes.
2315
2316
2317  */
2318
2319
2320  PyArrayObject *neighbours, *neighbour_edges,
2321    *normals, *edgelengths, *radii, *areas,
2322    *tri_full_flag,
2323    *stage_edge_values,
2324    *xmom_edge_values,
2325    *ymom_edge_values,
2326    *bed_edge_values,
2327    *stage_boundary_values,
2328    *xmom_boundary_values,
2329    *ymom_boundary_values,
2330    *stage_explicit_update,
2331    *xmom_explicit_update,
2332    *ymom_explicit_update,
2333    *already_computed_flux, //Tracks whether the flux across an edge has already been computed
2334    *max_speed_array; //Keeps track of max speeds for each triangle
2335
2336   
2337  double timestep, epsilon, H0, g;
2338  int optimise_dry_cells;
2339   
2340  // Convert Python arguments to C
2341  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOi",
2342            &timestep,
2343            &epsilon,
2344            &H0,
2345            &g,
2346            &neighbours,
2347            &neighbour_edges,
2348            &normals,
2349            &edgelengths, &radii, &areas,
2350            &tri_full_flag,
2351            &stage_edge_values,
2352            &xmom_edge_values,
2353            &ymom_edge_values,
2354            &bed_edge_values,
2355            &stage_boundary_values,
2356            &xmom_boundary_values,
2357            &ymom_boundary_values,
2358            &stage_explicit_update,
2359            &xmom_explicit_update,
2360            &ymom_explicit_update,
2361            &already_computed_flux,
2362            &max_speed_array,
2363            &optimise_dry_cells)) {
2364    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
2365    return NULL;
2366  }
2367
2368  // check that numpy array objects arrays are C contiguous memory
2369  CHECK_C_CONTIG(neighbours);
2370  CHECK_C_CONTIG(neighbour_edges);
2371  CHECK_C_CONTIG(normals);
2372  CHECK_C_CONTIG(edgelengths);
2373  CHECK_C_CONTIG(radii);
2374  CHECK_C_CONTIG(areas);
2375  CHECK_C_CONTIG(tri_full_flag);
2376  CHECK_C_CONTIG(stage_edge_values);
2377  CHECK_C_CONTIG(xmom_edge_values);
2378  CHECK_C_CONTIG(ymom_edge_values);
2379  CHECK_C_CONTIG(bed_edge_values);
2380  CHECK_C_CONTIG(stage_boundary_values);
2381  CHECK_C_CONTIG(xmom_boundary_values);
2382  CHECK_C_CONTIG(ymom_boundary_values);
2383  CHECK_C_CONTIG(stage_explicit_update);
2384  CHECK_C_CONTIG(xmom_explicit_update);
2385  CHECK_C_CONTIG(ymom_explicit_update);
2386  CHECK_C_CONTIG(already_computed_flux);
2387  CHECK_C_CONTIG(max_speed_array);
2388
2389  int number_of_elements = stage_edge_values -> dimensions[0];
2390
2391  // Call underlying flux computation routine and update
2392  // the explicit update arrays
2393  timestep = _compute_fluxes_central(number_of_elements,
2394                     timestep,
2395                     epsilon,
2396                     H0,
2397                     g,
2398                     (long*) neighbours -> data,
2399                     (long*) neighbour_edges -> data,
2400                     (double*) normals -> data,
2401                     (double*) edgelengths -> data, 
2402                     (double*) radii -> data, 
2403                     (double*) areas -> data,
2404                     (long*) tri_full_flag -> data,
2405                     (double*) stage_edge_values -> data,
2406                     (double*) xmom_edge_values -> data,
2407                     (double*) ymom_edge_values -> data,
2408                     (double*) bed_edge_values -> data,
2409                     (double*) stage_boundary_values -> data,
2410                     (double*) xmom_boundary_values -> data,
2411                     (double*) ymom_boundary_values -> data,
2412                     (double*) stage_explicit_update -> data,
2413                     (double*) xmom_explicit_update -> data,
2414                     (double*) ymom_explicit_update -> data,
2415                     (long*) already_computed_flux -> data,
2416                     (double*) max_speed_array -> data,
2417                     optimise_dry_cells);
2418 
2419  // Return updated flux timestep
2420  return Py_BuildValue("d", timestep);
2421}
2422
2423
2424
2425
2426
2427PyObject *compute_fluxes_ext_kinetic(PyObject *self, PyObject *args) {
2428  /*Compute all fluxes and the timestep suitable for all volumes
2429    in domain.
2430
2431    THIS IS AN EXPERIMENTAL FUNCTION - NORMALLY flux_function_central IS USED.
2432   
2433    Compute total flux for each conserved quantity using "flux_function_kinetic"
2434
2435    Fluxes across each edge are scaled by edgelengths and summed up
2436    Resulting flux is then scaled by area and stored in
2437    explicit_update for each of the three conserved quantities
2438    stage, xmomentum and ymomentum
2439
2440    The maximal allowable speed computed by the flux_function for each volume
2441    is converted to a timestep that must not be exceeded. The minimum of
2442    those is computed as the next overall timestep.
2443
2444    Python call:
2445    domain.timestep = compute_fluxes(timestep,
2446                                     domain.epsilon,
2447                                     domain.H0,
2448                                     domain.g,
2449                                     domain.neighbours,
2450                                     domain.neighbour_edges,
2451                                     domain.normals,
2452                                     domain.edgelengths,
2453                                     domain.radii,
2454                                     domain.areas,
2455                                     Stage.edge_values,
2456                                     Xmom.edge_values,
2457                                     Ymom.edge_values,
2458                                     Bed.edge_values,
2459                                     Stage.boundary_values,
2460                                     Xmom.boundary_values,
2461                                     Ymom.boundary_values,
2462                                     Stage.explicit_update,
2463                                     Xmom.explicit_update,
2464                                     Ymom.explicit_update,
2465                                     already_computed_flux)
2466
2467
2468    Post conditions:
2469      domain.explicit_update is reset to computed flux values
2470      domain.timestep is set to the largest step satisfying all volumes.
2471
2472
2473  */
2474
2475
2476  PyArrayObject *neighbours, *neighbour_edges,
2477    *normals, *edgelengths, *radii, *areas,
2478    *tri_full_flag,
2479    *stage_edge_values,
2480    *xmom_edge_values,
2481    *ymom_edge_values,
2482    *bed_edge_values,
2483    *stage_boundary_values,
2484    *xmom_boundary_values,
2485    *ymom_boundary_values,
2486    *stage_explicit_update,
2487    *xmom_explicit_update,
2488    *ymom_explicit_update,
2489    *already_computed_flux; // Tracks whether the flux across an edge has already been computed
2490
2491
2492  // Local variables
2493  double timestep, max_speed, epsilon, g, H0;
2494  double normal[2], ql[3], qr[3], zl, zr;
2495  double edgeflux[3]; //Work arrays for summing up fluxes
2496
2497  int number_of_elements, k, i, m, n;
2498  int ki, nm=0, ki2; //Index shorthands
2499  static long call=1;
2500
2501
2502  // Convert Python arguments to C
2503  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOO",
2504            &timestep,
2505            &epsilon,
2506            &H0,
2507            &g,
2508            &neighbours,
2509            &neighbour_edges,
2510            &normals,
2511            &edgelengths, &radii, &areas,
2512            &tri_full_flag,
2513            &stage_edge_values,
2514            &xmom_edge_values,
2515            &ymom_edge_values,
2516            &bed_edge_values,
2517            &stage_boundary_values,
2518            &xmom_boundary_values,
2519            &ymom_boundary_values,
2520            &stage_explicit_update,
2521            &xmom_explicit_update,
2522            &ymom_explicit_update,
2523            &already_computed_flux)) {
2524    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
2525    return NULL;
2526  }
2527  number_of_elements = stage_edge_values -> dimensions[0];
2528  call++;//a static local variable to which already_computed_flux is compared
2529  //set explicit_update to zero for all conserved_quantities.
2530  //This assumes compute_fluxes called before forcing terms
2531  for (k=0; k<number_of_elements; k++) {
2532    ((double *) stage_explicit_update -> data)[k]=0.0;
2533    ((double *) xmom_explicit_update -> data)[k]=0.0;
2534    ((double *) ymom_explicit_update -> data)[k]=0.0;
2535  }
2536  //Loop through neighbours and compute edge flux for each
2537  for (k=0; k<number_of_elements; k++) {
2538    for (i=0; i<3; i++) {
2539      ki = k*3+i;
2540      if (((long *) already_computed_flux->data)[ki]==call)//we've already computed the flux across this edge
2541    continue;
2542      ql[0] = ((double *) stage_edge_values -> data)[ki];
2543      ql[1] = ((double *) xmom_edge_values -> data)[ki];
2544      ql[2] = ((double *) ymom_edge_values -> data)[ki];
2545      zl =    ((double *) bed_edge_values -> data)[ki];
2546
2547      //Quantities at neighbour on nearest face
2548      n = ((long *) neighbours -> data)[ki];
2549      if (n < 0) {
2550    m = -n-1; //Convert negative flag to index
2551    qr[0] = ((double *) stage_boundary_values -> data)[m];
2552    qr[1] = ((double *) xmom_boundary_values -> data)[m];
2553    qr[2] = ((double *) ymom_boundary_values -> data)[m];
2554    zr = zl; //Extend bed elevation to boundary
2555      } else {
2556    m = ((long *) neighbour_edges -> data)[ki];
2557    nm = n*3+m;
2558    qr[0] = ((double *) stage_edge_values -> data)[nm];
2559    qr[1] = ((double *) xmom_edge_values -> data)[nm];
2560    qr[2] = ((double *) ymom_edge_values -> data)[nm];
2561    zr =    ((double *) bed_edge_values -> data)[nm];
2562      }
2563      // Outward pointing normal vector
2564      // normal = domain.normals[k, 2*i:2*i+2]
2565      ki2 = 2*ki; //k*6 + i*2
2566      normal[0] = ((double *) normals -> data)[ki2];
2567      normal[1] = ((double *) normals -> data)[ki2+1];
2568      //Edge flux computation
2569      flux_function_kinetic(ql, qr, zl, zr,
2570            normal[0], normal[1],
2571            epsilon, H0, g,
2572            edgeflux, &max_speed);
2573      //update triangle k
2574      ((long *) already_computed_flux->data)[ki]=call;
2575      ((double *) stage_explicit_update -> data)[k] -= edgeflux[0]*((double *) edgelengths -> data)[ki];
2576      ((double *) xmom_explicit_update -> data)[k] -= edgeflux[1]*((double *) edgelengths -> data)[ki];
2577      ((double *) ymom_explicit_update -> data)[k] -= edgeflux[2]*((double *) edgelengths -> data)[ki];
2578      //update the neighbour n
2579      if (n>=0){
2580    ((long *) already_computed_flux->data)[nm]=call;
2581    ((double *) stage_explicit_update -> data)[n] += edgeflux[0]*((double *) edgelengths -> data)[nm];
2582    ((double *) xmom_explicit_update -> data)[n] += edgeflux[1]*((double *) edgelengths -> data)[nm];
2583    ((double *) ymom_explicit_update -> data)[n] += edgeflux[2]*((double *) edgelengths -> data)[nm];
2584      }
2585      ///for (j=0; j<3; j++) {
2586    ///flux[j] -= edgeflux[j]*((double *) edgelengths -> data)[ki];
2587    ///}
2588    //Update timestep
2589    //timestep = min(timestep, domain.radii[k]/max_speed)
2590    //FIXME: SR Add parameter for CFL condition
2591    if ( ((long *) tri_full_flag -> data)[k] == 1) {
2592        if (max_speed > epsilon) {
2593            timestep = min(timestep, ((double *) radii -> data)[k]/max_speed);
2594            //maxspeed in flux_function is calculated as max(|u+a|,|u-a|)
2595            if (n>=0)
2596                timestep = min(timestep, ((double *) radii -> data)[n]/max_speed);
2597        }
2598    }
2599    } // end for i
2600    //Normalise by area and store for when all conserved
2601    //quantities get updated
2602    ((double *) stage_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
2603    ((double *) xmom_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
2604    ((double *) ymom_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
2605  } //end for k
2606  return Py_BuildValue("d", timestep);
2607}
2608
2609PyObject *protect(PyObject *self, PyObject *args) {
2610  //
2611  //    protect(minimum_allowed_height, maximum_allowed_speed, wc, zc, xmomc, ymomc)
2612
2613
2614  PyArrayObject
2615  *wc,            //Stage at centroids
2616  *zc,            //Elevation at centroids
2617  *xmomc,         //Momentums at centroids
2618  *ymomc;
2619
2620
2621  int N;
2622  double minimum_allowed_height, maximum_allowed_speed, epsilon;
2623
2624  // Convert Python arguments to C
2625  if (!PyArg_ParseTuple(args, "dddOOOO",
2626            &minimum_allowed_height,
2627            &maximum_allowed_speed,
2628            &epsilon,
2629            &wc, &zc, &xmomc, &ymomc)) {
2630    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: protect could not parse input arguments");
2631    return NULL;
2632  } 
2633
2634  N = wc -> dimensions[0];
2635
2636  _protect(N,
2637       minimum_allowed_height,
2638       maximum_allowed_speed,
2639       epsilon,
2640       (double*) wc -> data,
2641       (double*) zc -> data,
2642       (double*) xmomc -> data,
2643       (double*) ymomc -> data);
2644
2645  return Py_BuildValue("");
2646}
2647
2648
2649
2650PyObject *balance_deep_and_shallow(PyObject *self, PyObject *args) {
2651  // Compute linear combination between stage as computed by
2652  // gradient-limiters limiting using w, and stage computed by
2653  // gradient-limiters limiting using h (h-limiter).
2654  // The former takes precedence when heights are large compared to the
2655  // bed slope while the latter takes precedence when heights are
2656  // relatively small.  Anything in between is computed as a balanced
2657  // linear combination in order to avoid numerical disturbances which
2658  // would otherwise appear as a result of hard switching between
2659  // modes.
2660  //
2661  //    balance_deep_and_shallow(wc, zc, wv, zv,
2662  //                             xmomc, ymomc, xmomv, ymomv)
2663
2664
2665  PyArrayObject
2666    *wc,            //Stage at centroids
2667    *zc,            //Elevation at centroids
2668    *wv,            //Stage at vertices
2669    *zv,            //Elevation at vertices
2670    *hvbar,         //h-Limited depths at vertices
2671    *xmomc,         //Momentums at centroids and vertices
2672    *ymomc,
2673    *xmomv,
2674    *ymomv;
2675 
2676  PyObject *domain, *Tmp;
2677   
2678  double alpha_balance = 2.0;
2679  double H0;
2680
2681  int N, tight_slope_limiters, use_centroid_velocities; //, err;
2682
2683  // Convert Python arguments to C
2684  if (!PyArg_ParseTuple(args, "OOOOOOOOOO",
2685            &domain,
2686            &wc, &zc, 
2687            &wv, &zv, &hvbar,
2688            &xmomc, &ymomc, &xmomv, &ymomv)) {
2689    PyErr_SetString(PyExc_RuntimeError, 
2690            "shallow_water_ext.c: balance_deep_and_shallow could not parse input arguments");
2691    return NULL;
2692  } 
2693     
2694     
2695  // FIXME (Ole): I tested this without GetAttrString and got time down
2696  // marginally from 4.0s to 3.8s. Consider passing everything in
2697  // through ParseTuple and profile.
2698 
2699  // Pull out parameters
2700  Tmp = PyObject_GetAttrString(domain, "alpha_balance");
2701  if (!Tmp) {
2702    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not obtain object alpha_balance from domain");
2703    return NULL;
2704  } 
2705  alpha_balance = PyFloat_AsDouble(Tmp);
2706  Py_DECREF(Tmp);
2707
2708 
2709  Tmp = PyObject_GetAttrString(domain, "H0");
2710  if (!Tmp) {
2711    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not obtain object H0 from domain");
2712    return NULL;
2713  } 
2714  H0 = PyFloat_AsDouble(Tmp);
2715  Py_DECREF(Tmp);
2716
2717 
2718  Tmp = PyObject_GetAttrString(domain, "tight_slope_limiters");
2719  if (!Tmp) {
2720    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not obtain object tight_slope_limiters from domain");
2721    return NULL;
2722  } 
2723  tight_slope_limiters = PyInt_AsLong(Tmp);
2724  Py_DECREF(Tmp);
2725 
2726 
2727  Tmp = PyObject_GetAttrString(domain, "use_centroid_velocities");
2728  if (!Tmp) {
2729    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not obtain object use_centroid_velocities from domain");
2730    return NULL;
2731  } 
2732  use_centroid_velocities = PyInt_AsLong(Tmp);
2733  Py_DECREF(Tmp);
2734 
2735 
2736 
2737  N = wc -> dimensions[0];
2738  _balance_deep_and_shallow(N,
2739                (double*) wc -> data,
2740                (double*) zc -> data,
2741                (double*) wv -> data,
2742                (double*) zv -> data,
2743                (double*) hvbar -> data,
2744                (double*) xmomc -> data,
2745                (double*) ymomc -> data,
2746                (double*) xmomv -> data,
2747                (double*) ymomv -> data,
2748                H0,
2749                (int) tight_slope_limiters,
2750                (int) use_centroid_velocities,             
2751                alpha_balance);
2752
2753
2754  return Py_BuildValue("");
2755}
2756
2757
2758
2759
2760PyObject *assign_windfield_values(PyObject *self, PyObject *args) {
2761  //
2762  //      assign_windfield_values(xmom_update, ymom_update,
2763  //                              s_vec, phi_vec, self.const)
2764
2765
2766
2767  PyArrayObject   //(one element per triangle)
2768  *s_vec,         //Speeds
2769  *phi_vec,       //Bearings
2770  *xmom_update,   //Momentum updates
2771  *ymom_update;
2772
2773
2774  int N;
2775  double cw;
2776
2777  // Convert Python arguments to C
2778  if (!PyArg_ParseTuple(args, "OOOOd",
2779            &xmom_update,
2780            &ymom_update,
2781            &s_vec, &phi_vec,
2782            &cw)) {
2783    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: assign_windfield_values could not parse input arguments");
2784    return NULL;
2785  } 
2786           
2787
2788  N = xmom_update -> dimensions[0];
2789
2790  _assign_wind_field_values(N,
2791       (double*) xmom_update -> data,
2792       (double*) ymom_update -> data,
2793       (double*) s_vec -> data,
2794       (double*) phi_vec -> data,
2795       cw);
2796
2797  return Py_BuildValue("");
2798}
2799
2800
2801
2802
2803//-------------------------------
2804// Method table for python module
2805//-------------------------------
2806static struct PyMethodDef MethodTable[] = {
2807  /* The cast of the function is necessary since PyCFunction values
2808   * only take two PyObject* parameters, and rotate() takes
2809   * three.
2810   */
2811
2812  {"rotate", (PyCFunction)rotate, METH_VARARGS | METH_KEYWORDS, "Print out"},
2813  {"extrapolate_second_order_sw", extrapolate_second_order_sw, METH_VARARGS, "Print out"},
2814  {"compute_fluxes_ext_central", compute_fluxes_ext_central, METH_VARARGS, "Print out"},
2815  {"compute_fluxes_ext_central_new", compute_fluxes_ext_central_new, METH_VARARGS, "Print out"},
2816  {"compute_fluxes_ext_kinetic", compute_fluxes_ext_kinetic, METH_VARARGS, "Print out"},
2817  {"gravity", gravity, METH_VARARGS, "Print out"},
2818  {"manning_friction_old", manning_friction_old, METH_VARARGS, "Print out"},
2819  {"manning_friction_new", manning_friction_new, METH_VARARGS, "Print out"},
2820  {"flux_function_central", flux_function_central, METH_VARARGS, "Print out"}, 
2821  {"balance_deep_and_shallow", balance_deep_and_shallow,
2822   METH_VARARGS, "Print out"},
2823  {"protect", protect, METH_VARARGS | METH_KEYWORDS, "Print out"},
2824  {"assign_windfield_values", assign_windfield_values,
2825   METH_VARARGS | METH_KEYWORDS, "Print out"},
2826  {NULL, NULL}
2827};
2828
2829// Module initialisation
2830void initshallow_water_ext(void){
2831  Py_InitModule("shallow_water_ext", MethodTable);
2832
2833  import_array(); // Necessary for handling of NumPY structures
2834}
Note: See TracBrowser for help on using the repository browser.