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

Last change on this file since 4730 was 4730, checked in by ole, 17 years ago

Cosmetics

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
1262    on each triangle
1263   
1264    These values are calculated as follows:
1265    1) For each triangle not adjacent to a boundary, we consider the
1266       auxiliary triangle formed by the centroids of its three
1267       neighbours.
1268    2) For each conserved quantity, we integrate around the auxiliary
1269       triangle's boundary the product of the quantity and the outward
1270       normal vector. Dividing by the triangle area gives (a,b), the
1271       average of the vector (q_x,q_y) on the auxiliary triangle.
1272       We suppose that the linear reconstruction on the original
1273       triangle has gradient (a,b).
1274    3) Provisional vertex jumps dqv[0,1,2] are computed and these are
1275       then limited by calling the functions find_qmin_and_qmax and
1276       limit_gradient
1277
1278    Python call:
1279    extrapolate_second_order_sw(domain.surrogate_neighbours,
1280                                domain.number_of_boundaries
1281                                domain.centroid_coordinates,
1282                                Stage.centroid_values
1283                                Xmom.centroid_values
1284                                Ymom.centroid_values
1285                                domain.vertex_coordinates,
1286                                Stage.vertex_values,
1287                                Xmom.vertex_values,
1288                                Ymom.vertex_values)
1289
1290    Post conditions:
1291            The vertices of each triangle have values from a
1292            limited linear reconstruction
1293            based on centroid values
1294
1295  */
1296  PyArrayObject *surrogate_neighbours,
1297    *number_of_boundaries,
1298    *centroid_coordinates,
1299    *stage_centroid_values,
1300    *xmom_centroid_values,
1301    *ymom_centroid_values,
1302        *elevation_centroid_values,
1303    *vertex_coordinates,
1304    *stage_vertex_values,
1305    *xmom_vertex_values,
1306    *ymom_vertex_values,
1307        *elevation_vertex_values;
1308  PyObject *domain, *Tmp;
1309
1310 
1311  double beta_w, beta_w_dry, beta_uh, beta_uh_dry, beta_vh, beta_vh_dry;   
1312  double minimum_allowed_height, epsilon;
1313  int optimise_dry_cells, number_of_elements, e;
1314 
1315  // Provisional jumps from centroids to v'tices and safety factor re limiting
1316  // by which these jumps are limited
1317  // Convert Python arguments to C
1318  if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOi",
1319                        &domain,
1320                        &surrogate_neighbours,
1321                        &number_of_boundaries,
1322                        &centroid_coordinates,
1323                        &stage_centroid_values,
1324                        &xmom_centroid_values,
1325                        &ymom_centroid_values,
1326                        &elevation_centroid_values,
1327                        &vertex_coordinates,
1328                        &stage_vertex_values,
1329                        &xmom_vertex_values,
1330                        &ymom_vertex_values,
1331                        &elevation_vertex_values,
1332                        &optimise_dry_cells)) {                 
1333                       
1334    PyErr_SetString(PyExc_RuntimeError, 
1335                    "Input arguments to extrapolate_second_order_sw failed");
1336    return NULL;
1337  }
1338
1339  // Get the safety factor beta_w, set in the config.py file.
1340  // This is used in the limiting process
1341 
1342 
1343  Tmp = PyObject_GetAttrString(domain, "beta_w");
1344  if (!Tmp) {
1345    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_w from domain");
1346    return NULL;
1347  } 
1348  beta_w = PyFloat_AsDouble(Tmp);
1349  Py_DECREF(Tmp);
1350 
1351  Tmp = PyObject_GetAttrString(domain, "beta_w_dry");
1352  if (!Tmp) {
1353    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_w_dry from domain");
1354    return NULL;
1355  } 
1356  beta_w_dry = PyFloat_AsDouble(Tmp);
1357  Py_DECREF(Tmp);
1358 
1359  Tmp = PyObject_GetAttrString(domain, "beta_uh");
1360  if (!Tmp) {
1361    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_uh from domain");
1362    return NULL;
1363  } 
1364  beta_uh = PyFloat_AsDouble(Tmp);
1365  Py_DECREF(Tmp);
1366 
1367  Tmp = PyObject_GetAttrString(domain, "beta_uh_dry");
1368  if (!Tmp) {
1369    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_uh_dry from domain");
1370    return NULL;
1371  } 
1372  beta_uh_dry = PyFloat_AsDouble(Tmp);
1373  Py_DECREF(Tmp); 
1374
1375  Tmp = PyObject_GetAttrString(domain, "beta_vh");
1376  if (!Tmp) {
1377    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_vh from domain");
1378    return NULL;
1379  } 
1380  beta_vh = PyFloat_AsDouble(Tmp);
1381  Py_DECREF(Tmp);
1382 
1383  Tmp = PyObject_GetAttrString(domain, "beta_vh_dry");
1384  if (!Tmp) {
1385    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_vh_dry from domain");
1386    return NULL;
1387  } 
1388  beta_vh_dry = PyFloat_AsDouble(Tmp);
1389  Py_DECREF(Tmp);
1390 
1391  Tmp = PyObject_GetAttrString(domain, "minimum_allowed_height");
1392  if (!Tmp) {
1393    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object minimum_allowed_height");
1394    return NULL;
1395  } 
1396  minimum_allowed_height = PyFloat_AsDouble(Tmp);
1397  Py_DECREF(Tmp); 
1398
1399  Tmp = PyObject_GetAttrString(domain, "epsilon");
1400  if (!Tmp) {
1401    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object epsilon");
1402    return NULL;
1403  } 
1404  epsilon = PyFloat_AsDouble(Tmp);
1405  Py_DECREF(Tmp); 
1406 
1407 
1408  // Call underlying computational routine
1409  number_of_elements = stage_centroid_values -> dimensions[0]; 
1410  e = _extrapolate_second_order_sw(number_of_elements,
1411                                   epsilon,
1412                                   minimum_allowed_height,
1413                                   beta_w,
1414                                   beta_w_dry,
1415                                   beta_uh,
1416                                   beta_uh_dry,
1417                                   beta_vh,
1418                                   beta_vh_dry,
1419                                   (long*) surrogate_neighbours -> data,
1420                                   (long*) number_of_boundaries -> data,
1421                                   (double*) centroid_coordinates -> data,
1422                                   (double*) stage_centroid_values -> data,
1423                                   (double*) xmom_centroid_values -> data,
1424                                   (double*) ymom_centroid_values -> data,
1425                                   (double*) elevation_centroid_values -> data,
1426                                   (double*) vertex_coordinates -> data,
1427                                   (double*) stage_vertex_values -> data,
1428                                   (double*) xmom_vertex_values -> data,
1429                                   (double*) ymom_vertex_values -> data,
1430                                   (double*) elevation_vertex_values -> data,
1431                                   optimise_dry_cells);
1432  if (e == -1) {
1433    // Use error string set inside computational routine
1434    return NULL;
1435  }                               
1436 
1437 
1438  return Py_BuildValue("");
1439 
1440}// extrapolate_second-order_sw
1441
1442
1443
1444// FIXME (Ole): This function is obsolete as of 12 Sep 2007
1445PyObject *extrapolate_second_order_sw_original(PyObject *self, PyObject *args) {
1446  /*Compute the vertex values based on a linear reconstruction on each triangle
1447    These values are calculated as follows:
1448    1) For each triangle not adjacent to a boundary, we consider the auxiliary triangle
1449    formed by the centroids of its three neighbours.
1450    2) For each conserved quantity, we integrate around the auxiliary triangle's boundary the product
1451    of the quantity and the outward normal vector. Dividing by the triangle area gives (a,b), the average
1452    of the vector (q_x,q_y) on the auxiliary triangle. We suppose that the linear reconstruction on the
1453    original triangle has gradient (a,b).
1454    3) Provisional vertex jumps dqv[0,1,2] are computed and these are then limited by calling the functions
1455    find_qmin_and_qmax and limit_gradient
1456
1457    Python call:
1458    extrapolate_second_order_sw(domain.surrogate_neighbours,
1459                                domain.number_of_boundaries
1460                                domain.centroid_coordinates,
1461                                Stage.centroid_values
1462                                Xmom.centroid_values
1463                                Ymom.centroid_values
1464                                domain.vertex_coordinates,
1465                                Stage.vertex_values,
1466                                Xmom.vertex_values,
1467                                Ymom.vertex_values)
1468
1469    Post conditions:
1470            The vertices of each triangle have values from a limited linear reconstruction
1471            based on centroid values
1472
1473  */
1474  PyArrayObject *surrogate_neighbours,
1475    *number_of_boundaries,
1476    *centroid_coordinates,
1477    *stage_centroid_values,
1478    *xmom_centroid_values,
1479    *ymom_centroid_values,
1480        *elevation_centroid_values,
1481    *vertex_coordinates,
1482    *stage_vertex_values,
1483    *xmom_vertex_values,
1484    *ymom_vertex_values,
1485        *elevation_vertex_values;
1486  PyObject *domain, *Tmp;
1487  //Local variables
1488  double a, b;//gradient vector, not stored but used to calculate vertex values from centroids
1489  int number_of_elements,k,k0,k1,k2,k3,k6,coord_index,i;
1490  double x,y,x0,y0,x1,y1,x2,y2,xv0,yv0,xv1,yv1,xv2,yv2;//vertices of the auxiliary triangle
1491  double dx1,dx2,dy1,dy2,dxv0,dxv1,dxv2,dyv0,dyv1,dyv2,dq0,dq1,dq2,area2;
1492  double dqv[3], qmin, qmax, hmin, hmax;
1493  double hc, h0, h1, h2;
1494  double epsilon=1.0e-12;
1495  int optimise_dry_cells=1; // Optimisation flag   
1496  double beta_w, beta_w_dry, beta_uh, beta_uh_dry, beta_vh, beta_vh_dry, beta_tmp;
1497  double minimum_allowed_height;
1498 
1499  // Provisional jumps from centroids to v'tices and safety factor re limiting
1500  // by which these jumps are limited
1501  // Convert Python arguments to C
1502  if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOi",
1503                        &domain,
1504                        &surrogate_neighbours,
1505                        &number_of_boundaries,
1506                        &centroid_coordinates,
1507                        &stage_centroid_values,
1508                        &xmom_centroid_values,
1509                        &ymom_centroid_values,
1510                        &elevation_centroid_values,
1511                        &vertex_coordinates,
1512                        &stage_vertex_values,
1513                        &xmom_vertex_values,
1514                        &ymom_vertex_values,
1515                        &elevation_vertex_values,
1516                        &optimise_dry_cells)) {                 
1517                       
1518    PyErr_SetString(PyExc_RuntimeError, "Input arguments to extrapolate_second_order_sw failed");
1519    return NULL;
1520  }
1521
1522 
1523  // Get the safety factor beta_w, set in the config.py file. This is used in the limiting process
1524  Tmp = PyObject_GetAttrString(domain, "beta_w");
1525  if (!Tmp) {
1526    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_w from domain");
1527    return NULL;
1528  } 
1529  beta_w = PyFloat_AsDouble(Tmp);
1530  Py_DECREF(Tmp);
1531 
1532  Tmp = PyObject_GetAttrString(domain, "beta_w_dry");
1533  if (!Tmp) {
1534    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_w_dry from domain");
1535    return NULL;
1536  } 
1537  beta_w_dry = PyFloat_AsDouble(Tmp);
1538  Py_DECREF(Tmp);
1539 
1540  Tmp = PyObject_GetAttrString(domain, "beta_uh");
1541  if (!Tmp) {
1542    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_uh from domain");
1543    return NULL;
1544  } 
1545  beta_uh = PyFloat_AsDouble(Tmp);
1546  Py_DECREF(Tmp);
1547 
1548  Tmp = PyObject_GetAttrString(domain, "beta_uh_dry");
1549  if (!Tmp) {
1550    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_uh_dry from domain");
1551    return NULL;
1552  } 
1553  beta_uh_dry = PyFloat_AsDouble(Tmp);
1554  Py_DECREF(Tmp); 
1555
1556  Tmp = PyObject_GetAttrString(domain, "beta_vh");
1557  if (!Tmp) {
1558    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_vh from domain");
1559    return NULL;
1560  } 
1561  beta_vh = PyFloat_AsDouble(Tmp);
1562  Py_DECREF(Tmp);
1563 
1564  Tmp = PyObject_GetAttrString(domain, "beta_vh_dry");
1565  if (!Tmp) {
1566    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_vh_dry from domain");
1567    return NULL;
1568  } 
1569  beta_vh_dry = PyFloat_AsDouble(Tmp);
1570  Py_DECREF(Tmp);
1571 
1572  Tmp = PyObject_GetAttrString(domain, "minimum_allowed_height");
1573  if (!Tmp) {
1574    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object minimum_allowed_height");
1575    return NULL;
1576  } 
1577  minimum_allowed_height = PyFloat_AsDouble(Tmp);
1578  Py_DECREF(Tmp); 
1579
1580  Tmp = PyObject_GetAttrString(domain, "epsilon");
1581  if (!Tmp) {
1582    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object epsilon");
1583    return NULL;
1584  } 
1585  epsilon = PyFloat_AsDouble(Tmp);
1586  Py_DECREF(Tmp); 
1587   
1588  number_of_elements = stage_centroid_values -> dimensions[0];
1589  for (k=0; k<number_of_elements; k++) {
1590    k3=k*3;
1591    k6=k*6;
1592
1593   
1594    if (((long *) number_of_boundaries->data)[k]==3){
1595      // No neighbours, set gradient on the triangle to zero*/
1596      ((double *) stage_vertex_values->data)[k3]=((double *)stage_centroid_values->data)[k];
1597      ((double *) stage_vertex_values->data)[k3+1]=((double *)stage_centroid_values->data)[k];
1598      ((double *) stage_vertex_values->data)[k3+2]=((double *)stage_centroid_values->data)[k];
1599      ((double *) xmom_vertex_values->data)[k3]=((double *)xmom_centroid_values->data)[k];
1600      ((double *) xmom_vertex_values->data)[k3+1]=((double *)xmom_centroid_values->data)[k];
1601      ((double *) xmom_vertex_values->data)[k3+2]=((double *)xmom_centroid_values->data)[k];
1602      ((double *) ymom_vertex_values->data)[k3]=((double *)ymom_centroid_values->data)[k];
1603      ((double *) ymom_vertex_values->data)[k3+1]=((double *)ymom_centroid_values->data)[k];
1604      ((double *) ymom_vertex_values->data)[k3+2]=((double *)ymom_centroid_values->data)[k];
1605      continue;
1606    }
1607    else {
1608      // Triangle k has one or more neighbours.
1609      // Get centroid and vertex coordinates of the triangle
1610     
1611      // Get the vertex coordinates
1612      xv0=((double *)vertex_coordinates->data)[k6]; yv0=((double *)vertex_coordinates->data)[k6+1];
1613      xv1=((double *)vertex_coordinates->data)[k6+2]; yv1=((double *)vertex_coordinates->data)[k6+3];
1614      xv2=((double *)vertex_coordinates->data)[k6+4]; yv2=((double *)vertex_coordinates->data)[k6+5];
1615     
1616      // Get the centroid coordinates
1617      coord_index=2*k;
1618      x=((double *)centroid_coordinates->data)[coord_index];
1619      y=((double *)centroid_coordinates->data)[coord_index+1];
1620     
1621      // Store x- and y- differentials for the vertices of the FV triangle relative to the centroid
1622      dxv0=xv0-x; dxv1=xv1-x; dxv2=xv2-x;
1623      dyv0=yv0-y; dyv1=yv1-y; dyv2=yv2-y;
1624    }
1625
1626
1627   
1628   
1629   
1630           
1631    if (((long *)number_of_boundaries->data)[k]<=1){
1632   
1633      //==============================================
1634      // Number of boundaries <= 1
1635      //==============================================   
1636   
1637   
1638      // If no boundaries, auxiliary triangle is formed from the centroids of the three neighbours
1639      // If one boundary, auxiliary triangle is formed from this centroid and its two neighbours
1640      k0=((long *)surrogate_neighbours->data)[k3];
1641      k1=((long *)surrogate_neighbours->data)[k3+1];
1642      k2=((long *)surrogate_neighbours->data)[k3+2];
1643     
1644      // Get the auxiliary triangle's vertex coordinates (really the centroids of neighbouring triangles)
1645      coord_index=2*k0;
1646      x0=((double *)centroid_coordinates->data)[coord_index];
1647      y0=((double *)centroid_coordinates->data)[coord_index+1];
1648      coord_index=2*k1;
1649      x1=((double *)centroid_coordinates->data)[coord_index];
1650      y1=((double *)centroid_coordinates->data)[coord_index+1];
1651      coord_index=2*k2;
1652      x2=((double *)centroid_coordinates->data)[coord_index];
1653      y2=((double *)centroid_coordinates->data)[coord_index+1];
1654     
1655      // Store x- and y- differentials for the vertices of the auxiliary triangle
1656      dx1=x1-x0; dx2=x2-x0;
1657      dy1=y1-y0; dy2=y2-y0;
1658     
1659      // Calculate 2*area of the auxiliary triangle
1660      area2 = dy2*dx1 - dy1*dx2;//the triangle is guaranteed to be counter-clockwise
1661     
1662      // If the mesh is 'weird' near the boundary, the triangle might be flat or clockwise:
1663      if (area2<=0) {
1664        PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: negative triangle area encountered");
1665        return NULL;
1666      } 
1667     
1668      // Calculate heights of neighbouring cells
1669      hc = ((double *)stage_centroid_values->data)[k]  - ((double *)elevation_centroid_values->data)[k];
1670      h0 = ((double *)stage_centroid_values->data)[k0] - ((double *)elevation_centroid_values->data)[k0];
1671      h1 = ((double *)stage_centroid_values->data)[k1] - ((double *)elevation_centroid_values->data)[k1];
1672      h2 = ((double *)stage_centroid_values->data)[k2] - ((double *)elevation_centroid_values->data)[k2];
1673      hmin = min(hc,min(h0,min(h1,h2)));  // FIXME Don't need to include hc
1674     
1675     
1676      if (optimise_dry_cells) {     
1677        // Check if linear reconstruction is necessary for triangle k
1678        // This check will exclude dry cells.
1679
1680        hmax = max(h0,max(h1,h2));     
1681        if (hmax < epsilon) {
1682          continue;
1683        }
1684      }
1685
1686           
1687      //-----------------------------------
1688      // stage
1689      //-----------------------------------     
1690     
1691      // Calculate the difference between vertex 0 of the auxiliary triangle and the FV triangle centroid
1692      dq0=((double *)stage_centroid_values->data)[k0]-((double *)stage_centroid_values->data)[k];
1693     
1694      // Calculate differentials between the vertices of the auxiliary triangle
1695      dq1=((double *)stage_centroid_values->data)[k1]-((double *)stage_centroid_values->data)[k0];
1696      dq2=((double *)stage_centroid_values->data)[k2]-((double *)stage_centroid_values->data)[k0];
1697     
1698      // Calculate the gradient of stage on the auxiliary triangle
1699      a = dy2*dq1 - dy1*dq2;
1700      a /= area2;
1701      b = dx1*dq2 - dx2*dq1;
1702      b /= area2;
1703     
1704      // Calculate provisional jumps in stage from the centroid of the FV tri to its vertices, to be limited
1705      dqv[0]=a*dxv0+b*dyv0;
1706      dqv[1]=a*dxv1+b*dyv1;
1707      dqv[2]=a*dxv2+b*dyv2;
1708     
1709      // Now we want to find min and max of the centroid and the vertices of the auxiliary triangle
1710      // and compute jumps from the centroid to the min and max
1711      find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
1712     
1713      // Playing with dry wet interface
1714      hmin = qmin;
1715      beta_tmp = beta_w;
1716      if (hmin<minimum_allowed_height)
1717        beta_tmp = beta_w_dry;
1718      limit_gradient(dqv,qmin,qmax,beta_tmp);//the gradient will be limited
1719      for (i=0;i<3;i++)
1720        ((double *)stage_vertex_values->data)[k3+i]=((double *)stage_centroid_values->data)[k]+dqv[i];
1721     
1722     
1723      //-----------------------------------
1724      // xmomentum
1725      //-----------------------------------           
1726
1727      // Calculate the difference between vertex 0 of the auxiliary triangle and the FV triangle centroid
1728      dq0=((double *)xmom_centroid_values->data)[k0]-((double *)xmom_centroid_values->data)[k];
1729     
1730      // Calculate differentials between the vertices of the auxiliary triangle
1731      dq1=((double *)xmom_centroid_values->data)[k1]-((double *)xmom_centroid_values->data)[k0];
1732      dq2=((double *)xmom_centroid_values->data)[k2]-((double *)xmom_centroid_values->data)[k0];
1733     
1734      // Calculate the gradient of xmom on the auxiliary triangle
1735      a = dy2*dq1 - dy1*dq2;
1736      a /= area2;
1737      b = dx1*dq2 - dx2*dq1;
1738      b /= area2;
1739     
1740      // Calculate provisional jumps in stage from the centroid of the FV tri to its vertices, to be limited
1741      dqv[0]=a*dxv0+b*dyv0;
1742      dqv[1]=a*dxv1+b*dyv1;
1743      dqv[2]=a*dxv2+b*dyv2;
1744     
1745      // Now we want to find min and max of the centroid and the vertices of the auxiliary triangle
1746      // and compute jumps from the centroid to the min and max
1747      find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
1748      beta_tmp = beta_uh;
1749      if (hmin<minimum_allowed_height)
1750        beta_tmp = beta_uh_dry;
1751      limit_gradient(dqv,qmin,qmax,beta_tmp);//the gradient will be limited
1752      for (i=0;i<3;i++)
1753        ((double *)xmom_vertex_values->data)[k3+i]=((double *)xmom_centroid_values->data)[k]+dqv[i];
1754     
1755     
1756      //-----------------------------------
1757      // ymomentum
1758      //-----------------------------------                 
1759
1760      // Calculate the difference between vertex 0 of the auxiliary triangle and the FV triangle centroid
1761      dq0=((double *)ymom_centroid_values->data)[k0]-((double *)ymom_centroid_values->data)[k];
1762     
1763      // Calculate differentials between the vertices of the auxiliary triangle
1764      dq1=((double *)ymom_centroid_values->data)[k1]-((double *)ymom_centroid_values->data)[k0];
1765      dq2=((double *)ymom_centroid_values->data)[k2]-((double *)ymom_centroid_values->data)[k0];
1766     
1767      // Calculate the gradient of xmom on the auxiliary triangle
1768      a = dy2*dq1 - dy1*dq2;
1769      a /= area2;
1770      b = dx1*dq2 - dx2*dq1;
1771      b /= area2;
1772     
1773      // Calculate provisional jumps in stage from the centroid of the FV tri to its vertices, to be limited
1774      dqv[0]=a*dxv0+b*dyv0;
1775      dqv[1]=a*dxv1+b*dyv1;
1776      dqv[2]=a*dxv2+b*dyv2;
1777     
1778      // Now we want to find min and max of the centroid and the vertices of the auxiliary triangle
1779      // and compute jumps from the centroid to the min and max
1780      find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
1781      beta_tmp = beta_vh;
1782      if (hmin<minimum_allowed_height)
1783        beta_tmp = beta_vh_dry;
1784      limit_gradient(dqv,qmin,qmax,beta_tmp);//the gradient will be limited
1785      for (i=0;i<3;i++)
1786        ((double *)ymom_vertex_values->data)[k3+i]=((double *)ymom_centroid_values->data)[k]+dqv[i];
1787    } // End number_of_boundaries <=1
1788    else{
1789
1790      //==============================================
1791      // Number of boundaries == 2
1792      //==============================================       
1793       
1794      // One internal neighbour and gradient is in direction of the neighbour's centroid
1795     
1796      // Find the only internal neighbour
1797      for (k2=k3;k2<k3+3;k2++){//k2 just indexes the edges of triangle k
1798        if (((long *)surrogate_neighbours->data)[k2]!=k)//find internal neighbour of triabngle k
1799          break;
1800      }
1801      if ((k2==k3+3)) {//if we didn't find an internal neighbour
1802        PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: Internal neighbour not found");     
1803        return NULL;//error
1804      }
1805     
1806      k1=((long *)surrogate_neighbours->data)[k2];
1807     
1808      // The coordinates of the triangle are already (x,y). Get centroid of the neighbour (x1,y1)
1809      coord_index=2*k1;
1810      x1=((double *)centroid_coordinates->data)[coord_index];
1811      y1=((double *)centroid_coordinates->data)[coord_index+1];
1812     
1813      // Compute x- and y- distances between the centroid of the FV triangle and that of its neighbour
1814      dx1=x1-x; dy1=y1-y;
1815     
1816      // Set area2 as the square of the distance
1817      area2=dx1*dx1+dy1*dy1;
1818     
1819      // Set dx2=(x1-x0)/((x1-x0)^2+(y1-y0)^2) and dy2=(y1-y0)/((x1-x0)^2+(y1-y0)^2) which
1820      // respectively correspond to the x- and y- gradients of the conserved quantities
1821      dx2=1.0/area2;
1822      dy2=dx2*dy1;
1823      dx2*=dx1;
1824     
1825     
1826     
1827      //-----------------------------------
1828      // stage
1829      //-----------------------------------           
1830
1831      // Compute differentials
1832      dq1=((double *)stage_centroid_values->data)[k1]-((double *)stage_centroid_values->data)[k];
1833     
1834      // Calculate the gradient between the centroid of the FV triangle and that of its neighbour
1835      a=dq1*dx2;
1836      b=dq1*dy2;
1837     
1838      // Calculate provisional vertex jumps, to be limited
1839      dqv[0]=a*dxv0+b*dyv0;
1840      dqv[1]=a*dxv1+b*dyv1;
1841      dqv[2]=a*dxv2+b*dyv2;
1842     
1843      // Now limit the jumps
1844      if (dq1>=0.0){
1845        qmin=0.0;
1846        qmax=dq1;
1847      }
1848      else{
1849        qmin=dq1;
1850        qmax=0.0;
1851      }
1852     
1853     
1854      limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
1855      for (i=0;i<3;i++)
1856        ((double *)stage_vertex_values->data)[k3+i]=((double *)stage_centroid_values->data)[k]+dqv[i];
1857     
1858      //-----------------------------------
1859      // xmomentum
1860      //-----------------------------------                       
1861     
1862      // Compute differentials
1863      dq1=((double *)xmom_centroid_values->data)[k1]-((double *)xmom_centroid_values->data)[k];
1864     
1865      // Calculate the gradient between the centroid of the FV triangle and that of its neighbour
1866      a=dq1*dx2;
1867      b=dq1*dy2;
1868     
1869      // Calculate provisional vertex jumps, to be limited
1870      dqv[0]=a*dxv0+b*dyv0;
1871      dqv[1]=a*dxv1+b*dyv1;
1872      dqv[2]=a*dxv2+b*dyv2;
1873     
1874      // Now limit the jumps
1875      if (dq1>=0.0){
1876        qmin=0.0;
1877        qmax=dq1;
1878      }
1879      else{
1880        qmin=dq1;
1881        qmax=0.0;
1882      }
1883      limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
1884      for (i=0;i<3;i++)
1885        ((double *)xmom_vertex_values->data)[k3+i]=((double *)xmom_centroid_values->data)[k]+dqv[i];
1886     
1887      //-----------------------------------
1888      // ymomentum
1889      //-----------------------------------                       
1890
1891      // Compute differentials
1892      dq1=((double *)ymom_centroid_values->data)[k1]-((double *)ymom_centroid_values->data)[k];
1893     
1894      // Calculate the gradient between the centroid of the FV triangle and that of its neighbour
1895      a=dq1*dx2;
1896      b=dq1*dy2;
1897     
1898      // Calculate provisional vertex jumps, to be limited
1899      dqv[0]=a*dxv0+b*dyv0;
1900      dqv[1]=a*dxv1+b*dyv1;
1901      dqv[2]=a*dxv2+b*dyv2;
1902     
1903      // Now limit the jumps
1904      if (dq1>=0.0){
1905        qmin=0.0;
1906        qmax=dq1;
1907      }
1908      else{
1909        qmin=dq1;
1910        qmax=0.0;
1911      }
1912      limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
1913      for (i=0;i<3;i++)
1914        ((double *)ymom_vertex_values->data)[k3+i]=((double *)ymom_centroid_values->data)[k]+dqv[i];
1915    }//else [number_of_boundaries==2]
1916  }//for k=0 to number_of_elements-1
1917 
1918  return Py_BuildValue("");
1919}//extrapolate_second-order_sw
1920
1921
1922PyObject *rotate(PyObject *self, PyObject *args, PyObject *kwargs) {
1923  //
1924  // r = rotate(q, normal, direction=1)
1925  //
1926  // Where q is assumed to be a Float numeric array of length 3 and
1927  // normal a Float numeric array of length 2.
1928
1929  // FIXME(Ole): I don't think this is used anymore
1930
1931  PyObject *Q, *Normal;
1932  PyArrayObject *q, *r, *normal;
1933
1934  static char *argnames[] = {"q", "normal", "direction", NULL};
1935  int dimensions[1], i, direction=1;
1936  double n1, n2;
1937
1938  // Convert Python arguments to C
1939  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OO|i", argnames,
1940                                   &Q, &Normal, &direction)) {
1941    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: rotate could not parse input arguments");
1942    return NULL;
1943  } 
1944
1945  //Input checks (convert sequences into numeric arrays)
1946  q = (PyArrayObject *)
1947    PyArray_ContiguousFromObject(Q, PyArray_DOUBLE, 0, 0);
1948  normal = (PyArrayObject *)
1949    PyArray_ContiguousFromObject(Normal, PyArray_DOUBLE, 0, 0);
1950
1951
1952  if (normal -> dimensions[0] != 2) {
1953    PyErr_SetString(PyExc_RuntimeError, "Normal vector must have 2 components");
1954    return NULL;
1955  }
1956
1957  //Allocate space for return vector r (don't DECREF)
1958  dimensions[0] = 3;
1959  r = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
1960
1961  //Copy
1962  for (i=0; i<3; i++) {
1963    ((double *) (r -> data))[i] = ((double *) (q -> data))[i];
1964  }
1965
1966  //Get normal and direction
1967  n1 = ((double *) normal -> data)[0];
1968  n2 = ((double *) normal -> data)[1];
1969  if (direction == -1) n2 = -n2;
1970
1971  //Rotate
1972  _rotate((double *) r -> data, n1, n2);
1973
1974  //Release numeric arrays
1975  Py_DECREF(q);
1976  Py_DECREF(normal);
1977
1978  //return result using PyArray to avoid memory leak
1979  return PyArray_Return(r);
1980}
1981
1982
1983// Computational function for flux computation
1984double _compute_fluxes_central(int number_of_elements,
1985                               double timestep,
1986                               double epsilon,
1987                               double H0,
1988                               double g,
1989                               long* neighbours,
1990                               long* neighbour_edges,
1991                               double* normals,
1992                               double* edgelengths, 
1993                               double* radii, 
1994                               double* areas,
1995                               long* tri_full_flag,
1996                               double* stage_edge_values,
1997                               double* xmom_edge_values,
1998                               double* ymom_edge_values,
1999                               double* bed_edge_values,
2000                               double* stage_boundary_values,
2001                               double* xmom_boundary_values,
2002                               double* ymom_boundary_values,
2003                               double* stage_explicit_update,
2004                               double* xmom_explicit_update,
2005                               double* ymom_explicit_update,
2006                               long* already_computed_flux,
2007                               double* max_speed_array,
2008                               int optimise_dry_cells) {
2009                               
2010  // Local variables
2011  double max_speed, length, area, zl, zr;
2012  int k, i, m, n;
2013  int ki, nm=0, ki2; // Index shorthands
2014 
2015  // Workspace (making them static actually made function slightly slower (Ole)) 
2016  double ql[3], qr[3], edgeflux[3]; // Work array for summing up fluxes
2017
2018  static long call=1; // Static local variable flagging already computed flux
2019                               
2020       
2021  // Start computation
2022  call++; // Flag 'id' of flux calculation for this timestep
2023 
2024  // Set explicit_update to zero for all conserved_quantities.
2025  // This assumes compute_fluxes called before forcing terms
2026  for (k=0; k<number_of_elements; k++) {
2027    stage_explicit_update[k]=0.0;
2028    xmom_explicit_update[k]=0.0;
2029    ymom_explicit_update[k]=0.0; 
2030  }
2031
2032  // For all triangles
2033  for (k=0; k<number_of_elements; k++) {
2034   
2035    // Loop through neighbours and compute edge flux for each 
2036    for (i=0; i<3; i++) {
2037      ki = k*3+i; // Linear index (triangle k, edge i)
2038     
2039      if (already_computed_flux[ki] == call)
2040        // We've already computed the flux across this edge
2041        continue;
2042       
2043       
2044      ql[0] = stage_edge_values[ki];
2045      ql[1] = xmom_edge_values[ki];
2046      ql[2] = ymom_edge_values[ki];
2047      zl = bed_edge_values[ki];
2048
2049      // Quantities at neighbour on nearest face
2050      n = neighbours[ki];
2051      if (n < 0) {
2052        m = -n-1; // Convert negative flag to index
2053       
2054        qr[0] = stage_boundary_values[m];
2055        qr[1] = xmom_boundary_values[m];
2056        qr[2] = ymom_boundary_values[m];
2057        zr = zl; // Extend bed elevation to boundary
2058      } else {
2059        m = neighbour_edges[ki];
2060        nm = n*3+m; // Linear index (triangle n, edge m)
2061       
2062        qr[0] = stage_edge_values[nm];
2063        qr[1] = xmom_edge_values[nm];
2064        qr[2] = ymom_edge_values[nm];
2065        zr = bed_edge_values[nm];
2066      }
2067     
2068     
2069      if (optimise_dry_cells) {     
2070        // Check if flux calculation is necessary across this edge
2071        // This check will exclude dry cells.
2072        // This will also optimise cases where zl != zr as
2073        // long as both are dry
2074
2075        if ( fabs(ql[0] - zl) < epsilon && 
2076             fabs(qr[0] - zr) < epsilon ) {
2077          // Cell boundary is dry
2078         
2079          already_computed_flux[ki] = call; // #k Done 
2080          if (n>=0)
2081            already_computed_flux[nm] = call; // #n Done
2082       
2083          max_speed = 0.0;
2084          continue;
2085        }
2086      }
2087   
2088     
2089      // Outward pointing normal vector (domain.normals[k, 2*i:2*i+2])
2090      ki2 = 2*ki; //k*6 + i*2
2091
2092      // Edge flux computation (triangle k, edge i)
2093      flux_function_central(ql, qr, zl, zr,
2094                            normals[ki2], normals[ki2+1],
2095                            epsilon, H0, g,
2096                            edgeflux, &max_speed);
2097     
2098     
2099      // Multiply edgeflux by edgelength
2100      length = edgelengths[ki];
2101      edgeflux[0] *= length;           
2102      edgeflux[1] *= length;           
2103      edgeflux[2] *= length;                       
2104     
2105     
2106      // Update triangle k with flux from edge i
2107      stage_explicit_update[k] -= edgeflux[0];
2108      xmom_explicit_update[k] -= edgeflux[1];
2109      ymom_explicit_update[k] -= edgeflux[2];
2110     
2111      already_computed_flux[ki] = call; // #k Done
2112     
2113     
2114      // Update neighbour n with same flux but reversed sign
2115      if (n>=0) {
2116        stage_explicit_update[n] += edgeflux[0];
2117        xmom_explicit_update[n] += edgeflux[1];
2118        ymom_explicit_update[n] += edgeflux[2];
2119       
2120        already_computed_flux[nm] = call; // #n Done
2121      }
2122
2123
2124      // Update timestep based on edge i and possibly neighbour n
2125      if (tri_full_flag[k] == 1) {
2126        if (max_speed > epsilon) {
2127          timestep = min(timestep, radii[k]/max_speed);
2128          if (n>=0) 
2129            timestep = min(timestep, radii[n]/max_speed);
2130        }
2131      }
2132     
2133    } // End edge i
2134   
2135   
2136    // Normalise triangle k by area and store for when all conserved
2137    // quantities get updated
2138    area = areas[k];
2139    stage_explicit_update[k] /= area;
2140    xmom_explicit_update[k] /= area;
2141    ymom_explicit_update[k] /= area;
2142   
2143   
2144    // Keep track of maximal speeds
2145    max_speed_array[k] = max_speed;   
2146   
2147  } // End triangle k
2148       
2149                               
2150                               
2151  return timestep;
2152}
2153
2154
2155PyObject *compute_fluxes_ext_central(PyObject *self, PyObject *args) {
2156  /*Compute all fluxes and the timestep suitable for all volumes
2157    in domain.
2158
2159    Compute total flux for each conserved quantity using "flux_function_central"
2160
2161    Fluxes across each edge are scaled by edgelengths and summed up
2162    Resulting flux is then scaled by area and stored in
2163    explicit_update for each of the three conserved quantities
2164    stage, xmomentum and ymomentum
2165
2166    The maximal allowable speed computed by the flux_function for each volume
2167    is converted to a timestep that must not be exceeded. The minimum of
2168    those is computed as the next overall timestep.
2169
2170    Python call:
2171    domain.timestep = compute_fluxes(timestep,
2172                                     domain.epsilon,
2173                                     domain.H0,
2174                                     domain.g,
2175                                     domain.neighbours,
2176                                     domain.neighbour_edges,
2177                                     domain.normals,
2178                                     domain.edgelengths,
2179                                     domain.radii,
2180                                     domain.areas,
2181                                     tri_full_flag,
2182                                     Stage.edge_values,
2183                                     Xmom.edge_values,
2184                                     Ymom.edge_values,
2185                                     Bed.edge_values,
2186                                     Stage.boundary_values,
2187                                     Xmom.boundary_values,
2188                                     Ymom.boundary_values,
2189                                     Stage.explicit_update,
2190                                     Xmom.explicit_update,
2191                                     Ymom.explicit_update,
2192                                     already_computed_flux,
2193                                     optimise_dry_cells)                                     
2194
2195
2196    Post conditions:
2197      domain.explicit_update is reset to computed flux values
2198      domain.timestep is set to the largest step satisfying all volumes.
2199
2200
2201  */
2202
2203
2204  PyArrayObject *neighbours, *neighbour_edges,
2205    *normals, *edgelengths, *radii, *areas,
2206    *tri_full_flag,
2207    *stage_edge_values,
2208    *xmom_edge_values,
2209    *ymom_edge_values,
2210    *bed_edge_values,
2211    *stage_boundary_values,
2212    *xmom_boundary_values,
2213    *ymom_boundary_values,
2214    *stage_explicit_update,
2215    *xmom_explicit_update,
2216    *ymom_explicit_update,
2217    *already_computed_flux, //Tracks whether the flux across an edge has already been computed
2218    *max_speed_array; //Keeps track of max speeds for each triangle
2219
2220   
2221  double timestep, epsilon, H0, g;
2222  int optimise_dry_cells;
2223   
2224  // Convert Python arguments to C
2225  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOi",
2226                        &timestep,
2227                        &epsilon,
2228                        &H0,
2229                        &g,
2230                        &neighbours,
2231                        &neighbour_edges,
2232                        &normals,
2233                        &edgelengths, &radii, &areas,
2234                        &tri_full_flag,
2235                        &stage_edge_values,
2236                        &xmom_edge_values,
2237                        &ymom_edge_values,
2238                        &bed_edge_values,
2239                        &stage_boundary_values,
2240                        &xmom_boundary_values,
2241                        &ymom_boundary_values,
2242                        &stage_explicit_update,
2243                        &xmom_explicit_update,
2244                        &ymom_explicit_update,
2245                        &already_computed_flux,
2246                        &max_speed_array,
2247                        &optimise_dry_cells)) {
2248    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
2249    return NULL;
2250  }
2251
2252 
2253  int number_of_elements = stage_edge_values -> dimensions[0];
2254
2255  // Call underlying flux computation routine and update
2256  // the explicit update arrays
2257  timestep = _compute_fluxes_central(number_of_elements,
2258                                     timestep,
2259                                     epsilon,
2260                                     H0,
2261                                     g,
2262                                     (long*) neighbours -> data,
2263                                     (long*) neighbour_edges -> data,
2264                                     (double*) normals -> data,
2265                                     (double*) edgelengths -> data, 
2266                                     (double*) radii -> data, 
2267                                     (double*) areas -> data,
2268                                     (long*) tri_full_flag -> data,
2269                                     (double*) stage_edge_values -> data,
2270                                     (double*) xmom_edge_values -> data,
2271                                     (double*) ymom_edge_values -> data,
2272                                     (double*) bed_edge_values -> data,
2273                                     (double*) stage_boundary_values -> data,
2274                                     (double*) xmom_boundary_values -> data,
2275                                     (double*) ymom_boundary_values -> data,
2276                                     (double*) stage_explicit_update -> data,
2277                                     (double*) xmom_explicit_update -> data,
2278                                     (double*) ymom_explicit_update -> data,
2279                                     (long*) already_computed_flux -> data,
2280                                     (double*) max_speed_array -> data,
2281                                     optimise_dry_cells);
2282 
2283  // Return updated flux timestep
2284  return Py_BuildValue("d", timestep);
2285}
2286
2287
2288
2289
2290// THIS FUNCTION IS NOW OBSOLETE
2291PyObject *compute_fluxes_ext_central_original(PyObject *self, PyObject *args) {
2292  /*Compute all fluxes and the timestep suitable for all volumes
2293    in domain.
2294
2295    Compute total flux for each conserved quantity using "flux_function_central"
2296
2297    Fluxes across each edge are scaled by edgelengths and summed up
2298    Resulting flux is then scaled by area and stored in
2299    explicit_update for each of the three conserved quantities
2300    stage, xmomentum and ymomentum
2301
2302    The maximal allowable speed computed by the flux_function for each volume
2303    is converted to a timestep that must not be exceeded. The minimum of
2304    those is computed as the next overall timestep.
2305
2306    Python call:
2307    domain.timestep = compute_fluxes(timestep,
2308                                     domain.epsilon,
2309                                     domain.H0,
2310                                     domain.g,
2311                                     domain.neighbours,
2312                                     domain.neighbour_edges,
2313                                     domain.normals,
2314                                     domain.edgelengths,
2315                                     domain.radii,
2316                                     domain.areas,
2317                                     tri_full_flag,
2318                                     Stage.edge_values,
2319                                     Xmom.edge_values,
2320                                     Ymom.edge_values,
2321                                     Bed.edge_values,
2322                                     Stage.boundary_values,
2323                                     Xmom.boundary_values,
2324                                     Ymom.boundary_values,
2325                                     Stage.explicit_update,
2326                                     Xmom.explicit_update,
2327                                     Ymom.explicit_update,
2328                                     already_computed_flux,
2329                                     optimise_dry_cells)                                     
2330
2331
2332    Post conditions:
2333      domain.explicit_update is reset to computed flux values
2334      domain.timestep is set to the largest step satisfying all volumes.
2335
2336
2337  */
2338
2339
2340  PyArrayObject *neighbours, *neighbour_edges,
2341    *normals, *edgelengths, *radii, *areas,
2342    *tri_full_flag,
2343    *stage_edge_values,
2344    *xmom_edge_values,
2345    *ymom_edge_values,
2346    *bed_edge_values,
2347    *stage_boundary_values,
2348    *xmom_boundary_values,
2349    *ymom_boundary_values,
2350    *stage_explicit_update,
2351    *xmom_explicit_update,
2352    *ymom_explicit_update,
2353    *already_computed_flux, //Tracks whether the flux across an edge has already been computed
2354    *max_speed_array; //Keeps track of max speeds for each triangle
2355
2356
2357  // Local variables
2358  double timestep, max_speed, epsilon, g, H0, length, area;
2359  int optimise_dry_cells=0; // Optimisation flag 
2360  double normal[2], ql[3], qr[3], zl, zr;
2361  double edgeflux[3]; // Work array for summing up fluxes
2362
2363  int number_of_elements, k, i, m, n;
2364
2365  int ki, nm=0, ki2; // Index shorthands
2366  static long call=1; // Static local variable flagging already computed flux
2367
2368
2369  // Convert Python arguments to C
2370  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOi",
2371                        &timestep,
2372                        &epsilon,
2373                        &H0,
2374                        &g,
2375                        &neighbours,
2376                        &neighbour_edges,
2377                        &normals,
2378                        &edgelengths, &radii, &areas,
2379                        &tri_full_flag,
2380                        &stage_edge_values,
2381                        &xmom_edge_values,
2382                        &ymom_edge_values,
2383                        &bed_edge_values,
2384                        &stage_boundary_values,
2385                        &xmom_boundary_values,
2386                        &ymom_boundary_values,
2387                        &stage_explicit_update,
2388                        &xmom_explicit_update,
2389                        &ymom_explicit_update,
2390                        &already_computed_flux,
2391                        &max_speed_array,
2392                        &optimise_dry_cells)) {
2393    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
2394    return NULL;
2395  }
2396 
2397 
2398  number_of_elements = stage_edge_values -> dimensions[0];
2399
2400  call++; // Flag 'id' of flux calculation for this timestep
2401 
2402  // Set explicit_update to zero for all conserved_quantities.
2403  // This assumes compute_fluxes called before forcing terms
2404  for (k=0; k<number_of_elements; k++) {
2405    ((double *) stage_explicit_update -> data)[k]=0.0;
2406    ((double *) xmom_explicit_update -> data)[k]=0.0;
2407    ((double *) ymom_explicit_update -> data)[k]=0.0; 
2408  }
2409 
2410  // For all triangles
2411  for (k=0; k<number_of_elements; k++) {
2412 
2413    // Loop through neighbours and compute edge flux for each 
2414    for (i=0; i<3; i++) {
2415      ki = k*3+i; // Linear index (triangle k, edge i)
2416     
2417      if (((long *) already_computed_flux->data)[ki] == call)
2418        // We've already computed the flux across this edge
2419        continue;
2420       
2421       
2422      ql[0] = ((double *) stage_edge_values -> data)[ki];
2423      ql[1] = ((double *) xmom_edge_values -> data)[ki];
2424      ql[2] = ((double *) ymom_edge_values -> data)[ki];
2425      zl =    ((double *) bed_edge_values -> data)[ki];
2426
2427      // Quantities at neighbour on nearest face
2428      n = ((long *) neighbours -> data)[ki];
2429      if (n < 0) {
2430        m = -n-1; // Convert negative flag to index
2431       
2432        qr[0] = ((double *) stage_boundary_values -> data)[m];
2433        qr[1] = ((double *) xmom_boundary_values -> data)[m];
2434        qr[2] = ((double *) ymom_boundary_values -> data)[m];
2435        zr = zl; //Extend bed elevation to boundary
2436      } else {
2437        m = ((long *) neighbour_edges -> data)[ki];
2438        nm = n*3+m; // Linear index (triangle n, edge m)
2439       
2440        qr[0] = ((double *) stage_edge_values -> data)[nm];
2441        qr[1] = ((double *) xmom_edge_values -> data)[nm];
2442        qr[2] = ((double *) ymom_edge_values -> data)[nm];
2443        zr =    ((double *) bed_edge_values -> data)[nm];
2444      }
2445     
2446     
2447      if (optimise_dry_cells) {     
2448        // Check if flux calculation is necessary across this edge
2449        // This check will exclude dry cells.
2450        // This will also optimise cases where zl != zr as
2451        // long as both are dry
2452
2453        if ( fabs(ql[0] - zl) < epsilon && 
2454             fabs(qr[0] - zr) < epsilon ) {
2455          // Cell boundary is dry
2456         
2457          ((long *) already_computed_flux -> data)[ki] = call; // #k Done       
2458          if (n>=0)
2459            ((long *) already_computed_flux -> data)[nm] = call; // #n Done
2460       
2461          max_speed = 0.0;
2462          continue;
2463        }
2464      }
2465     
2466           
2467      // Outward pointing normal vector (domain.normals[k, 2*i:2*i+2])
2468      ki2 = 2*ki; //k*6 + i*2
2469      normal[0] = ((double *) normals -> data)[ki2];
2470      normal[1] = ((double *) normals -> data)[ki2+1];
2471     
2472
2473      // Edge flux computation (triangle k, edge i)
2474      flux_function_central(ql, qr, zl, zr,
2475                            normal[0], normal[1],
2476                            epsilon, H0, g,
2477                            edgeflux, &max_speed);
2478     
2479     
2480      // Multiply edgeflux by edgelength
2481      length = ((double *) edgelengths -> data)[ki];
2482      edgeflux[0] *= length;           
2483      edgeflux[1] *= length;           
2484      edgeflux[2] *= length;                       
2485     
2486     
2487      // Update triangle k with flux from edge i
2488      ((double *) stage_explicit_update -> data)[k] -= edgeflux[0];
2489      ((double *) xmom_explicit_update -> data)[k] -= edgeflux[1];
2490      ((double *) ymom_explicit_update -> data)[k] -= edgeflux[2];
2491     
2492      ((long *) already_computed_flux -> data)[ki] = call; // #k Done
2493     
2494     
2495      // Update neighbour n with same flux but reversed sign
2496      if (n>=0){
2497        ((double *) stage_explicit_update -> data)[n] += edgeflux[0];
2498        ((double *) xmom_explicit_update -> data)[n] += edgeflux[1];
2499        ((double *) ymom_explicit_update -> data)[n] += edgeflux[2];
2500       
2501        ((long *) already_computed_flux -> data)[nm] = call; // #n Done
2502      }
2503     
2504     
2505      // Update timestep based on edge i and possibly neighbour n
2506      if ( ((long *) tri_full_flag -> data)[k] == 1) {
2507        if (max_speed > epsilon) {
2508          timestep = min(timestep, ((double *) radii -> data)[k]/max_speed);
2509          if (n>=0) 
2510            timestep = min(timestep, ((double *) radii -> data)[n]/max_speed);
2511        }
2512      }
2513     
2514    } // End edge i
2515   
2516   
2517    // Normalise triangle k by area and store for when all conserved
2518    // quantities get updated
2519    area = ((double *) areas -> data)[k];
2520    ((double *) stage_explicit_update -> data)[k] /= area;
2521    ((double *) xmom_explicit_update -> data)[k] /= area;
2522    ((double *) ymom_explicit_update -> data)[k] /= area;
2523   
2524   
2525    // Keep track of maximal speeds
2526    ((double *) max_speed_array -> data)[k] = max_speed;   
2527   
2528  } // End triangle k
2529 
2530  return Py_BuildValue("d", timestep);
2531}
2532
2533
2534PyObject *compute_fluxes_ext_kinetic(PyObject *self, PyObject *args) {
2535  /*Compute all fluxes and the timestep suitable for all volumes
2536    in domain.
2537
2538    Compute total flux for each conserved quantity using "flux_function_central"
2539
2540    Fluxes across each edge are scaled by edgelengths and summed up
2541    Resulting flux is then scaled by area and stored in
2542    explicit_update for each of the three conserved quantities
2543    stage, xmomentum and ymomentum
2544
2545    The maximal allowable speed computed by the flux_function for each volume
2546    is converted to a timestep that must not be exceeded. The minimum of
2547    those is computed as the next overall timestep.
2548
2549    Python call:
2550    domain.timestep = compute_fluxes(timestep,
2551                                     domain.epsilon,
2552                                     domain.H0,
2553                                     domain.g,
2554                                     domain.neighbours,
2555                                     domain.neighbour_edges,
2556                                     domain.normals,
2557                                     domain.edgelengths,
2558                                     domain.radii,
2559                                     domain.areas,
2560                                     Stage.edge_values,
2561                                     Xmom.edge_values,
2562                                     Ymom.edge_values,
2563                                     Bed.edge_values,
2564                                     Stage.boundary_values,
2565                                     Xmom.boundary_values,
2566                                     Ymom.boundary_values,
2567                                     Stage.explicit_update,
2568                                     Xmom.explicit_update,
2569                                     Ymom.explicit_update,
2570                                     already_computed_flux)
2571
2572
2573    Post conditions:
2574      domain.explicit_update is reset to computed flux values
2575      domain.timestep is set to the largest step satisfying all volumes.
2576
2577
2578  */
2579
2580
2581  PyArrayObject *neighbours, *neighbour_edges,
2582    *normals, *edgelengths, *radii, *areas,
2583    *tri_full_flag,
2584    *stage_edge_values,
2585    *xmom_edge_values,
2586    *ymom_edge_values,
2587    *bed_edge_values,
2588    *stage_boundary_values,
2589    *xmom_boundary_values,
2590    *ymom_boundary_values,
2591    *stage_explicit_update,
2592    *xmom_explicit_update,
2593    *ymom_explicit_update,
2594    *already_computed_flux;//tracks whether the flux across an edge has already been computed
2595
2596
2597  //Local variables
2598  double timestep, max_speed, epsilon, g, H0;
2599  double normal[2], ql[3], qr[3], zl, zr;
2600  double edgeflux[3]; //Work arrays for summing up fluxes
2601
2602  int number_of_elements, k, i, m, n;
2603  int ki, nm=0, ki2; //Index shorthands
2604  static long call=1;
2605
2606
2607  // Convert Python arguments to C
2608  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOO",
2609                        &timestep,
2610                        &epsilon,
2611                        &H0,
2612                        &g,
2613                        &neighbours,
2614                        &neighbour_edges,
2615                        &normals,
2616                        &edgelengths, &radii, &areas,
2617                        &tri_full_flag,
2618                        &stage_edge_values,
2619                        &xmom_edge_values,
2620                        &ymom_edge_values,
2621                        &bed_edge_values,
2622                        &stage_boundary_values,
2623                        &xmom_boundary_values,
2624                        &ymom_boundary_values,
2625                        &stage_explicit_update,
2626                        &xmom_explicit_update,
2627                        &ymom_explicit_update,
2628                        &already_computed_flux)) {
2629    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
2630    return NULL;
2631  }
2632  number_of_elements = stage_edge_values -> dimensions[0];
2633  call++;//a static local variable to which already_computed_flux is compared
2634  //set explicit_update to zero for all conserved_quantities.
2635  //This assumes compute_fluxes called before forcing terms
2636  for (k=0; k<number_of_elements; k++) {
2637    ((double *) stage_explicit_update -> data)[k]=0.0;
2638    ((double *) xmom_explicit_update -> data)[k]=0.0;
2639    ((double *) ymom_explicit_update -> data)[k]=0.0;
2640  }
2641  //Loop through neighbours and compute edge flux for each
2642  for (k=0; k<number_of_elements; k++) {
2643    for (i=0; i<3; i++) {
2644      ki = k*3+i;
2645      if (((long *) already_computed_flux->data)[ki]==call)//we've already computed the flux across this edge
2646        continue;
2647      ql[0] = ((double *) stage_edge_values -> data)[ki];
2648      ql[1] = ((double *) xmom_edge_values -> data)[ki];
2649      ql[2] = ((double *) ymom_edge_values -> data)[ki];
2650      zl =    ((double *) bed_edge_values -> data)[ki];
2651
2652      //Quantities at neighbour on nearest face
2653      n = ((long *) neighbours -> data)[ki];
2654      if (n < 0) {
2655        m = -n-1; //Convert negative flag to index
2656        qr[0] = ((double *) stage_boundary_values -> data)[m];
2657        qr[1] = ((double *) xmom_boundary_values -> data)[m];
2658        qr[2] = ((double *) ymom_boundary_values -> data)[m];
2659        zr = zl; //Extend bed elevation to boundary
2660      } else {
2661        m = ((long *) neighbour_edges -> data)[ki];
2662        nm = n*3+m;
2663        qr[0] = ((double *) stage_edge_values -> data)[nm];
2664        qr[1] = ((double *) xmom_edge_values -> data)[nm];
2665        qr[2] = ((double *) ymom_edge_values -> data)[nm];
2666        zr =    ((double *) bed_edge_values -> data)[nm];
2667      }
2668      // Outward pointing normal vector
2669      // normal = domain.normals[k, 2*i:2*i+2]
2670      ki2 = 2*ki; //k*6 + i*2
2671      normal[0] = ((double *) normals -> data)[ki2];
2672      normal[1] = ((double *) normals -> data)[ki2+1];
2673      //Edge flux computation
2674      flux_function_kinetic(ql, qr, zl, zr,
2675                    normal[0], normal[1],
2676                    epsilon, H0, g,
2677                    edgeflux, &max_speed);
2678      //update triangle k
2679      ((long *) already_computed_flux->data)[ki]=call;
2680      ((double *) stage_explicit_update -> data)[k] -= edgeflux[0]*((double *) edgelengths -> data)[ki];
2681      ((double *) xmom_explicit_update -> data)[k] -= edgeflux[1]*((double *) edgelengths -> data)[ki];
2682      ((double *) ymom_explicit_update -> data)[k] -= edgeflux[2]*((double *) edgelengths -> data)[ki];
2683      //update the neighbour n
2684      if (n>=0){
2685        ((long *) already_computed_flux->data)[nm]=call;
2686        ((double *) stage_explicit_update -> data)[n] += edgeflux[0]*((double *) edgelengths -> data)[nm];
2687        ((double *) xmom_explicit_update -> data)[n] += edgeflux[1]*((double *) edgelengths -> data)[nm];
2688        ((double *) ymom_explicit_update -> data)[n] += edgeflux[2]*((double *) edgelengths -> data)[nm];
2689      }
2690      ///for (j=0; j<3; j++) {
2691        ///flux[j] -= edgeflux[j]*((double *) edgelengths -> data)[ki];
2692        ///}
2693        //Update timestep
2694        //timestep = min(timestep, domain.radii[k]/max_speed)
2695        //FIXME: SR Add parameter for CFL condition
2696    if ( ((long *) tri_full_flag -> data)[k] == 1) {
2697            if (max_speed > epsilon) {
2698                timestep = min(timestep, ((double *) radii -> data)[k]/max_speed);
2699                //maxspeed in flux_function is calculated as max(|u+a|,|u-a|)
2700                if (n>=0)
2701                    timestep = min(timestep, ((double *) radii -> data)[n]/max_speed);
2702            }
2703    }
2704    } // end for i
2705    //Normalise by area and store for when all conserved
2706    //quantities get updated
2707    ((double *) stage_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
2708    ((double *) xmom_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
2709    ((double *) ymom_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
2710  } //end for k
2711  return Py_BuildValue("d", timestep);
2712}
2713
2714PyObject *protect(PyObject *self, PyObject *args) {
2715  //
2716  //    protect(minimum_allowed_height, maximum_allowed_speed, wc, zc, xmomc, ymomc)
2717
2718
2719  PyArrayObject
2720  *wc,            //Stage at centroids
2721  *zc,            //Elevation at centroids
2722  *xmomc,         //Momentums at centroids
2723  *ymomc;
2724
2725
2726  int N;
2727  double minimum_allowed_height, maximum_allowed_speed, epsilon;
2728
2729  // Convert Python arguments to C
2730  if (!PyArg_ParseTuple(args, "dddOOOO",
2731                        &minimum_allowed_height,
2732                        &maximum_allowed_speed,
2733                        &epsilon,
2734                        &wc, &zc, &xmomc, &ymomc)) {
2735    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: protect could not parse input arguments");
2736    return NULL;
2737  } 
2738
2739  N = wc -> dimensions[0];
2740
2741  _protect(N,
2742           minimum_allowed_height,
2743           maximum_allowed_speed,
2744           epsilon,
2745           (double*) wc -> data,
2746           (double*) zc -> data,
2747           (double*) xmomc -> data,
2748           (double*) ymomc -> data);
2749
2750  return Py_BuildValue("");
2751}
2752
2753
2754
2755PyObject *balance_deep_and_shallow(PyObject *self, PyObject *args) {
2756  //
2757  //    balance_deep_and_shallow(wc, zc, hc, wv, zv, hv,
2758  //                             xmomc, ymomc, xmomv, ymomv)
2759
2760
2761  PyArrayObject
2762    *wc,            //Stage at centroids
2763    *zc,            //Elevation at centroids
2764    *hc,            //Height at centroids
2765    *wv,            //Stage at vertices
2766    *zv,            //Elevation at vertices
2767    *hv,            //Depths at vertices
2768    *hvbar,         //h-Limited depths at vertices
2769    *xmomc,         //Momentums at centroids and vertices
2770    *ymomc,
2771    *xmomv,
2772    *ymomv;
2773 
2774  PyObject *domain, *Tmp;
2775   
2776  double alpha_balance = 2.0;
2777  double H0;
2778
2779  int N, tight_slope_limiters; //, err;
2780
2781  // Convert Python arguments to C
2782  if (!PyArg_ParseTuple(args, "OOOOOOOOOOOO",
2783                        &domain,
2784                        &wc, &zc, &hc,
2785                        &wv, &zv, &hv, &hvbar,
2786                        &xmomc, &ymomc, &xmomv, &ymomv)) {
2787    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not parse input arguments");
2788    return NULL;
2789  } 
2790         
2791         
2792  // FIXME (Ole): I tested this without GetAttrString and got time down
2793  // marginally from 4.0s to 3.8s. Consider passing everything in
2794  // through ParseTuple and profile.
2795 
2796  // Pull out parameters
2797  Tmp = PyObject_GetAttrString(domain, "alpha_balance");
2798  if (!Tmp) {
2799    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not obtain object alpha_balance from domain");
2800    return NULL;
2801  } 
2802  alpha_balance = PyFloat_AsDouble(Tmp);
2803  Py_DECREF(Tmp);
2804
2805 
2806  Tmp = PyObject_GetAttrString(domain, "H0");
2807  if (!Tmp) {
2808    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not obtain object H0 from domain");
2809    return NULL;
2810  } 
2811  H0 = PyFloat_AsDouble(Tmp);
2812  Py_DECREF(Tmp);
2813
2814 
2815  Tmp = PyObject_GetAttrString(domain, "tight_slope_limiters");
2816  if (!Tmp) {
2817    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not obtain object tight_slope_limiters from domain");
2818    return NULL;
2819  } 
2820  tight_slope_limiters = PyInt_AsLong(Tmp);
2821  Py_DECREF(Tmp);
2822 
2823 
2824     
2825  //alpha_balance = 2.0;
2826  //H0 = 0.001;
2827  //tight_slope_limiters = 1;
2828 
2829  N = wc -> dimensions[0];
2830
2831  _balance_deep_and_shallow(N,
2832                            (double*) wc -> data,
2833                            (double*) zc -> data,
2834                            (double*) hc -> data,
2835                            (double*) wv -> data,
2836                            (double*) zv -> data,
2837                            (double*) hv -> data,
2838                            (double*) hvbar -> data,
2839                            (double*) xmomc -> data,
2840                            (double*) ymomc -> data,
2841                            (double*) xmomv -> data,
2842                            (double*) ymomv -> data,
2843                            H0,
2844                            (int) tight_slope_limiters,
2845                            alpha_balance);
2846
2847
2848  return Py_BuildValue("");
2849}
2850
2851
2852
2853PyObject *h_limiter(PyObject *self, PyObject *args) {
2854
2855  PyObject *domain, *Tmp;
2856  PyArrayObject
2857    *hv, *hc, //Depth at vertices and centroids
2858    *hvbar,   //Limited depth at vertices (return values)
2859    *neighbours;
2860
2861  int k, i, n, N, k3;
2862  int dimensions[2];
2863  double beta_h; //Safety factor (see config.py)
2864  double *hmin, *hmax, hn;
2865
2866  // Convert Python arguments to C
2867  if (!PyArg_ParseTuple(args, "OOO", &domain, &hc, &hv)) {
2868    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: h_limiter could not parse input arguments");
2869    return NULL;
2870  } 
2871 
2872  neighbours = get_consecutive_array(domain, "neighbours");
2873
2874  //Get safety factor beta_h
2875  Tmp = PyObject_GetAttrString(domain, "beta_h");
2876  if (!Tmp) {
2877    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: h_limiter could not obtain object beta_h from domain");
2878    return NULL;
2879  } 
2880  beta_h = PyFloat_AsDouble(Tmp);
2881
2882  Py_DECREF(Tmp);
2883
2884  N = hc -> dimensions[0];
2885
2886  //Create hvbar
2887  dimensions[0] = N;
2888  dimensions[1] = 3;
2889  hvbar = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
2890
2891
2892  //Find min and max of this and neighbour's centroid values
2893  hmin = malloc(N * sizeof(double));
2894  hmax = malloc(N * sizeof(double));
2895  for (k=0; k<N; k++) {
2896    k3=k*3;
2897
2898    hmin[k] = ((double*) hc -> data)[k];
2899    hmax[k] = hmin[k];
2900
2901    for (i=0; i<3; i++) {
2902      n = ((long*) neighbours -> data)[k3+i];
2903
2904      //Initialise hvbar with values from hv
2905      ((double*) hvbar -> data)[k3+i] = ((double*) hv -> data)[k3+i];
2906
2907      if (n >= 0) {
2908        hn = ((double*) hc -> data)[n]; //Neighbour's centroid value
2909
2910        hmin[k] = min(hmin[k], hn);
2911        hmax[k] = max(hmax[k], hn);
2912      }
2913    }
2914  }
2915
2916  // Call underlying standard routine
2917  _limit(N, beta_h, (double*) hc -> data, (double*) hvbar -> data, hmin, hmax);
2918
2919  // // //Py_DECREF(domain); //FIXME: NEcessary?
2920  free(hmin);
2921  free(hmax);
2922
2923  //return result using PyArray to avoid memory leak
2924  return PyArray_Return(hvbar);
2925  //return Py_BuildValue("");
2926}
2927
2928PyObject *h_limiter_sw(PyObject *self, PyObject *args) {
2929  //a faster version of h_limiter above
2930  PyObject *domain, *Tmp;
2931  PyArrayObject
2932    *hv, *hc, //Depth at vertices and centroids
2933    *hvbar,   //Limited depth at vertices (return values)
2934    *neighbours;
2935
2936  int k, i, N, k3,k0,k1,k2;
2937  int dimensions[2];
2938  double beta_h; //Safety factor (see config.py)
2939  double hmin, hmax, dh[3];
2940// Convert Python arguments to C
2941  if (!PyArg_ParseTuple(args, "OOO", &domain, &hc, &hv)) {
2942    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: h_limiter_sw could not parse input arguments");
2943    return NULL;
2944  } 
2945  neighbours = get_consecutive_array(domain, "neighbours");
2946
2947  //Get safety factor beta_h
2948  Tmp = PyObject_GetAttrString(domain, "beta_h");
2949  if (!Tmp) {
2950    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: h_limiter_sw could not obtain object beta_h from domain");
2951    return NULL;
2952  } 
2953  beta_h = PyFloat_AsDouble(Tmp);
2954
2955  Py_DECREF(Tmp);
2956
2957  N = hc -> dimensions[0];
2958
2959  //Create hvbar
2960  dimensions[0] = N;
2961  dimensions[1] = 3;
2962  hvbar = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
2963  for (k=0;k<N;k++){
2964    k3=k*3;
2965    //get the ids of the neighbours
2966    k0 = ((long*) neighbours -> data)[k3];
2967    k1 = ((long*) neighbours -> data)[k3+1];
2968    k2 = ((long*) neighbours -> data)[k3+2];
2969    //set hvbar provisionally
2970    for (i=0;i<3;i++){
2971      ((double*) hvbar -> data)[k3+i] = ((double*) hv -> data)[k3+i];
2972      dh[i]=((double*) hvbar -> data)[k3+i]-((double*) hc -> data)[k];
2973    }
2974    hmin=((double*) hc -> data)[k];
2975    hmax=hmin;
2976    if (k0>=0){
2977      hmin=min(hmin,((double*) hc -> data)[k0]);
2978      hmax=max(hmax,((double*) hc -> data)[k0]);
2979    }
2980    if (k1>=0){
2981      hmin=min(hmin,((double*) hc -> data)[k1]);
2982      hmax=max(hmax,((double*) hc -> data)[k1]);
2983    }
2984    if (k2>=0){
2985      hmin=min(hmin,((double*) hc -> data)[k2]);
2986      hmax=max(hmax,((double*) hc -> data)[k2]);
2987    }
2988    hmin-=((double*) hc -> data)[k];
2989    hmax-=((double*) hc -> data)[k];
2990    limit_gradient(dh,hmin,hmax,beta_h);
2991    for (i=0;i<3;i++)
2992      ((double*) hvbar -> data)[k3+i] = ((double*) hc -> data)[k]+dh[i];
2993  }
2994  return PyArray_Return(hvbar);
2995}
2996
2997PyObject *assign_windfield_values(PyObject *self, PyObject *args) {
2998  //
2999  //      assign_windfield_values(xmom_update, ymom_update,
3000  //                              s_vec, phi_vec, self.const)
3001
3002
3003
3004  PyArrayObject   //(one element per triangle)
3005  *s_vec,         //Speeds
3006  *phi_vec,       //Bearings
3007  *xmom_update,   //Momentum updates
3008  *ymom_update;
3009
3010
3011  int N;
3012  double cw;
3013
3014  // Convert Python arguments to C
3015  if (!PyArg_ParseTuple(args, "OOOOd",
3016                        &xmom_update,
3017                        &ymom_update,
3018                        &s_vec, &phi_vec,
3019                        &cw)) {
3020    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: assign_windfield_values could not parse input arguments");
3021    return NULL;
3022  } 
3023                       
3024
3025  N = xmom_update -> dimensions[0];
3026
3027  _assign_wind_field_values(N,
3028           (double*) xmom_update -> data,
3029           (double*) ymom_update -> data,
3030           (double*) s_vec -> data,
3031           (double*) phi_vec -> data,
3032           cw);
3033
3034  return Py_BuildValue("");
3035}
3036
3037
3038
3039
3040//////////////////////////////////////////
3041// Method table for python module
3042static struct PyMethodDef MethodTable[] = {
3043  /* The cast of the function is necessary since PyCFunction values
3044   * only take two PyObject* parameters, and rotate() takes
3045   * three.
3046   */
3047
3048  {"rotate", (PyCFunction)rotate, METH_VARARGS | METH_KEYWORDS, "Print out"},
3049  {"extrapolate_second_order_sw", extrapolate_second_order_sw, METH_VARARGS, "Print out"},
3050  {"compute_fluxes_ext_central", compute_fluxes_ext_central, METH_VARARGS, "Print out"},
3051  {"compute_fluxes_ext_kinetic", compute_fluxes_ext_kinetic, METH_VARARGS, "Print out"},
3052  {"gravity", gravity, METH_VARARGS, "Print out"},
3053  {"manning_friction", manning_friction, METH_VARARGS, "Print out"},
3054  {"balance_deep_and_shallow", balance_deep_and_shallow,
3055   METH_VARARGS, "Print out"},
3056  {"h_limiter", h_limiter,
3057   METH_VARARGS, "Print out"},
3058  {"h_limiter_sw", h_limiter_sw,
3059   METH_VARARGS, "Print out"},
3060  {"protect", protect, METH_VARARGS | METH_KEYWORDS, "Print out"},
3061  {"assign_windfield_values", assign_windfield_values,
3062   METH_VARARGS | METH_KEYWORDS, "Print out"},
3063  //{"distribute_to_vertices_and_edges",
3064  // distribute_to_vertices_and_edges, METH_VARARGS},
3065  //{"update_conserved_quantities",
3066  // update_conserved_quantities, METH_VARARGS},
3067  //{"set_initialcondition",
3068  // set_initialcondition, METH_VARARGS},
3069  {NULL, NULL}
3070};
3071
3072// Module initialisation
3073void initshallow_water_ext(void){
3074  Py_InitModule("shallow_water_ext", MethodTable);
3075
3076  import_array();     //Necessary for handling of NumPY structures
3077}
Note: See TracBrowser for help on using the repository browser.