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

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

More minor optimisations of flux calculations.

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