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

Last change on this file since 5531 was 5442, checked in by ole, 16 years ago

Retired h-limiter and beta_h as per ticket:194.
All unit tests and validation tests pass.

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