source: anuga_core/source/anuga/shallow_water/shallow_water_ext.c @ 4729

Last change on this file since 4729 was 4729, checked in by ole, 16 years ago

Refactored the function extrapolate_second_order_sw into a gateway routine and a computational routine.

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