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

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

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