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

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

Separated flux function into Gateway and Computational part.
This makes the code more readable and also improved the speed from
8.60s to 7.97s in the okushiri performance example on nautilus.
That is an 8.6% improvement for the flux computation and a 2% overall
improvement.

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