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

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

More small optimisations

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