source: branches/numpy/anuga/shallow_water/shallow_water_ext.c @ 6451

Last change on this file since 6451 was 6410, checked in by rwilson, 15 years ago

numpy changes.

File size: 74.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 "numpy/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  // check that numpy array objects arrays are C contiguous memory
882  CHECK_C_CONTIG(h);
883  CHECK_C_CONTIG(v);
884  CHECK_C_CONTIG(x);
885  CHECK_C_CONTIG(xmom);
886  CHECK_C_CONTIG(ymom);
887
888  N = h -> dimensions[0];
889  for (k=0; k<N; k++) {
890    k3 = 3*k;  // base index
891
892    // Get bathymetry
893    z0 = ((double*) v -> data)[k3 + 0];
894    z1 = ((double*) v -> data)[k3 + 1];
895    z2 = ((double*) v -> data)[k3 + 2];
896
897    // Optimise for flat bed
898    // Note (Ole): This didn't produce measurable speed up.
899    // Revisit later
900    //if (fabs(z0-z1)<epsilon && fabs(z1-z2)<epsilon) {
901    //  continue;
902    //}
903
904    // Get average depth from centroid values
905    avg_h = ((double *) h -> data)[k];
906
907    // Compute bed slope
908    k6 = 6*k;  // base index
909 
910    x0 = ((double*) x -> data)[k6 + 0];
911    y0 = ((double*) x -> data)[k6 + 1];
912    x1 = ((double*) x -> data)[k6 + 2];
913    y1 = ((double*) x -> data)[k6 + 3];
914    x2 = ((double*) x -> data)[k6 + 4];
915    y2 = ((double*) x -> data)[k6 + 5];
916
917
918    _gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2, &zx, &zy);
919
920    // Update momentum
921    ((double*) xmom -> data)[k] += -g*zx*avg_h;
922    ((double*) ymom -> data)[k] += -g*zy*avg_h;
923  }
924
925  return Py_BuildValue("");
926}
927
928
929PyObject *manning_friction(PyObject *self, PyObject *args) {
930  //
931  // manning_friction(g, eps, h, uh, vh, eta, xmom_update, ymom_update)
932  //
933
934
935  PyArrayObject *w, *z, *uh, *vh, *eta, *xmom, *ymom;
936  int N;
937  double g, eps;
938
939  if (!PyArg_ParseTuple(args, "ddOOOOOOO",
940                        &g, &eps, &w, &z, &uh, &vh, &eta,
941                        &xmom, &ymom)) {
942    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: manning_friction could not parse input arguments");
943    return NULL;
944  }
945
946  // check that numpy array objects arrays are C contiguous memory
947  CHECK_C_CONTIG(w);
948  CHECK_C_CONTIG(z);
949  CHECK_C_CONTIG(uh);
950  CHECK_C_CONTIG(vh);
951  CHECK_C_CONTIG(eta);
952  CHECK_C_CONTIG(xmom);
953  CHECK_C_CONTIG(ymom);
954
955  N = w -> dimensions[0];
956  _manning_friction(g, eps, N,
957                    (double*) w -> data,
958                    (double*) z -> data,
959                    (double*) uh -> data,
960                    (double*) vh -> data,
961                    (double*) eta -> data,
962                    (double*) xmom -> data,
963                    (double*) ymom -> data);
964
965  return Py_BuildValue("");
966}
967
968
969/*
970PyObject *manning_friction_explicit(PyObject *self, PyObject *args) {
971  //
972  // manning_friction_explicit(g, eps, h, uh, vh, eta, xmom_update, ymom_update)
973  //
974
975
976  PyArrayObject *w, *z, *uh, *vh, *eta, *xmom, *ymom;
977  int N;
978  double g, eps;
979
980  if (!PyArg_ParseTuple(args, "ddOOOOOOO",
981                        &g, &eps, &w, &z, &uh, &vh, &eta,
982                        &xmom, &ymom))
983    return NULL;
984
985  // check that numpy array objects arrays are C contiguous memory
986  CHECK_C_CONTIG(w);
987  CHECK_C_CONTIG(z);
988  CHECK_C_CONTIG(uh);
989  CHECK_C_CONTIG(vh);
990  CHECK_C_CONTIG(eta);
991  CHECK_C_CONTIG(xmom);
992  CHECK_C_CONTIG(ymom);
993
994  N = w -> dimensions[0];
995  _manning_friction_explicit(g, eps, N,
996                    (double*) w -> data,
997                    (double*) z -> data,
998                    (double*) uh -> data,
999                    (double*) vh -> data,
1000                    (double*) eta -> data,
1001                    (double*) xmom -> data,
1002                    (double*) ymom -> data);
1003
1004  return Py_BuildValue("");
1005}
1006*/
1007
1008
1009
1010// Computational routine
1011int _extrapolate_second_order_sw(int number_of_elements,
1012                                  double epsilon,
1013                                  double minimum_allowed_height,
1014                                  double beta_w,
1015                                  double beta_w_dry,
1016                                  double beta_uh,
1017                                  double beta_uh_dry,
1018                                  double beta_vh,
1019                                  double beta_vh_dry,
1020                                  long* surrogate_neighbours,
1021                                  long* number_of_boundaries,
1022                                  double* centroid_coordinates,
1023                                  double* stage_centroid_values,
1024                                  double* xmom_centroid_values,
1025                                  double* ymom_centroid_values,
1026                                  double* elevation_centroid_values,
1027                                  double* vertex_coordinates,
1028                                  double* stage_vertex_values,
1029                                  double* xmom_vertex_values,
1030                                  double* ymom_vertex_values,
1031                                  double* elevation_vertex_values,
1032                                  int optimise_dry_cells) {
1033                                 
1034                                 
1035
1036  // Local variables
1037  double a, b; // Gradient vector used to calculate vertex values from centroids
1038  int k,k0,k1,k2,k3,k6,coord_index,i;
1039  double x,y,x0,y0,x1,y1,x2,y2,xv0,yv0,xv1,yv1,xv2,yv2; // Vertices of the auxiliary triangle
1040  double dx1,dx2,dy1,dy2,dxv0,dxv1,dxv2,dyv0,dyv1,dyv2,dq0,dq1,dq2,area2;
1041  double dqv[3], qmin, qmax, hmin, hmax;
1042  double hc, h0, h1, h2, beta_tmp, hfactor;
1043
1044       
1045  for (k=0; k<number_of_elements; k++) {
1046    k3=k*3;
1047    k6=k*6;
1048
1049   
1050    if (number_of_boundaries[k]==3){
1051      // No neighbours, set gradient on the triangle to zero
1052     
1053      stage_vertex_values[k3]   = stage_centroid_values[k];
1054      stage_vertex_values[k3+1] = stage_centroid_values[k];
1055      stage_vertex_values[k3+2] = stage_centroid_values[k];
1056      xmom_vertex_values[k3]    = xmom_centroid_values[k];
1057      xmom_vertex_values[k3+1]  = xmom_centroid_values[k];
1058      xmom_vertex_values[k3+2]  = xmom_centroid_values[k];
1059      ymom_vertex_values[k3]    = ymom_centroid_values[k];
1060      ymom_vertex_values[k3+1]  = ymom_centroid_values[k];
1061      ymom_vertex_values[k3+2]  = ymom_centroid_values[k];
1062     
1063      continue;
1064    }
1065    else {
1066      // Triangle k has one or more neighbours.
1067      // Get centroid and vertex coordinates of the triangle
1068     
1069      // Get the vertex coordinates
1070      xv0 = vertex_coordinates[k6];   yv0=vertex_coordinates[k6+1];
1071      xv1 = vertex_coordinates[k6+2]; yv1=vertex_coordinates[k6+3];
1072      xv2 = vertex_coordinates[k6+4]; yv2=vertex_coordinates[k6+5];
1073     
1074      // Get the centroid coordinates
1075      coord_index=2*k;
1076      x=centroid_coordinates[coord_index];
1077      y=centroid_coordinates[coord_index+1];
1078     
1079      // Store x- and y- differentials for the vertices of
1080      // triangle k relative to the centroid
1081      dxv0=xv0-x; dxv1=xv1-x; dxv2=xv2-x;
1082      dyv0=yv0-y; dyv1=yv1-y; dyv2=yv2-y;
1083    }
1084
1085
1086           
1087    if (number_of_boundaries[k]<=1){
1088   
1089      //==============================================
1090      // Number of boundaries <= 1
1091      //==============================================   
1092   
1093   
1094      // If no boundaries, auxiliary triangle is formed
1095      // from the centroids of the three neighbours
1096      // If one boundary, auxiliary triangle is formed
1097      // from this centroid and its two neighbours
1098     
1099      k0=surrogate_neighbours[k3];
1100      k1=surrogate_neighbours[k3+1];
1101      k2=surrogate_neighbours[k3+2];
1102     
1103      // Get the auxiliary triangle's vertex coordinates
1104      // (really the centroids of neighbouring triangles)
1105      coord_index=2*k0;
1106      x0=centroid_coordinates[coord_index];
1107      y0=centroid_coordinates[coord_index+1];
1108     
1109      coord_index=2*k1;
1110      x1=centroid_coordinates[coord_index];
1111      y1=centroid_coordinates[coord_index+1];
1112     
1113      coord_index=2*k2;
1114      x2=centroid_coordinates[coord_index];
1115      y2=centroid_coordinates[coord_index+1];
1116     
1117      // Store x- and y- differentials for the vertices
1118      // of the auxiliary triangle
1119      dx1=x1-x0; dx2=x2-x0;
1120      dy1=y1-y0; dy2=y2-y0;
1121     
1122      // Calculate 2*area of the auxiliary triangle
1123      // The triangle is guaranteed to be counter-clockwise     
1124      area2 = dy2*dx1 - dy1*dx2; 
1125     
1126      // If the mesh is 'weird' near the boundary,
1127      // the triangle might be flat or clockwise:
1128      if (area2<=0) {
1129        PyErr_SetString(PyExc_RuntimeError, 
1130                        "shallow_water_ext.c: negative triangle area encountered");
1131        return -1;
1132      } 
1133     
1134      // Calculate heights of neighbouring cells
1135      hc = stage_centroid_values[k]  - elevation_centroid_values[k];
1136      h0 = stage_centroid_values[k0] - elevation_centroid_values[k0];
1137      h1 = stage_centroid_values[k1] - elevation_centroid_values[k1];
1138      h2 = stage_centroid_values[k2] - elevation_centroid_values[k2];
1139      hmin = min(min(h0,min(h1,h2)),hc);
1140      //hfactor = hc/(hc + 1.0);
1141
1142      hfactor = 0.0;
1143      if (hmin > 0.001 ) {
1144          hfactor = (hmin-0.001)/(hmin+0.004);
1145      }
1146     
1147      if (optimise_dry_cells) {     
1148        // Check if linear reconstruction is necessary for triangle k
1149        // This check will exclude dry cells.
1150
1151        hmax = max(h0,max(h1,h2));     
1152        if (hmax < epsilon) {
1153          continue;
1154        }
1155      }
1156
1157           
1158      //-----------------------------------
1159      // stage
1160      //-----------------------------------     
1161     
1162      // Calculate the difference between vertex 0 of the auxiliary
1163      // triangle and the centroid of triangle k
1164      dq0=stage_centroid_values[k0]-stage_centroid_values[k];
1165     
1166      // Calculate differentials between the vertices
1167      // of the auxiliary triangle (centroids of neighbouring triangles)
1168      dq1=stage_centroid_values[k1]-stage_centroid_values[k0];
1169      dq2=stage_centroid_values[k2]-stage_centroid_values[k0];
1170     
1171      // Calculate the gradient of stage on the auxiliary triangle
1172      a = dy2*dq1 - dy1*dq2;
1173      a /= area2;
1174      b = dx1*dq2 - dx2*dq1;
1175      b /= area2;
1176     
1177      // Calculate provisional jumps in stage from the centroid
1178      // of triangle k to its vertices, to be limited
1179      dqv[0]=a*dxv0+b*dyv0;
1180      dqv[1]=a*dxv1+b*dyv1;
1181      dqv[2]=a*dxv2+b*dyv2;
1182     
1183      // Now we want to find min and max of the centroid and the
1184      // vertices of the auxiliary triangle and compute jumps
1185      // from the centroid to the min and max
1186      find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
1187     
1188      // Playing with dry wet interface
1189      //hmin = qmin;
1190      //beta_tmp = beta_w_dry;
1191      //if (hmin>minimum_allowed_height)
1192      beta_tmp = beta_w_dry + (beta_w - beta_w_dry) * hfactor;
1193       
1194      //printf("min_alled_height = %f\n",minimum_allowed_height);
1195      //printf("hmin = %f\n",hmin);
1196      //printf("beta_w = %f\n",beta_w);
1197      //printf("beta_tmp = %f\n",beta_tmp);
1198      // Limit the gradient
1199      limit_gradient(dqv,qmin,qmax,beta_tmp); 
1200     
1201      for (i=0;i<3;i++)
1202        stage_vertex_values[k3+i]=stage_centroid_values[k]+dqv[i];
1203     
1204     
1205      //-----------------------------------
1206      // xmomentum
1207      //-----------------------------------           
1208
1209      // Calculate the difference between vertex 0 of the auxiliary
1210      // triangle and the centroid of triangle k     
1211      dq0=xmom_centroid_values[k0]-xmom_centroid_values[k];
1212     
1213      // Calculate differentials between the vertices
1214      // of the auxiliary triangle
1215      dq1=xmom_centroid_values[k1]-xmom_centroid_values[k0];
1216      dq2=xmom_centroid_values[k2]-xmom_centroid_values[k0];
1217     
1218      // Calculate the gradient of xmom on the auxiliary triangle
1219      a = dy2*dq1 - dy1*dq2;
1220      a /= area2;
1221      b = dx1*dq2 - dx2*dq1;
1222      b /= area2;
1223     
1224      // Calculate provisional jumps in stage from the centroid
1225      // of triangle k to its vertices, to be limited     
1226      dqv[0]=a*dxv0+b*dyv0;
1227      dqv[1]=a*dxv1+b*dyv1;
1228      dqv[2]=a*dxv2+b*dyv2;
1229     
1230      // Now we want to find min and max of the centroid and the
1231      // vertices of the auxiliary triangle and compute jumps
1232      // from the centroid to the min and max
1233      find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
1234      //beta_tmp = beta_uh;
1235      //if (hmin<minimum_allowed_height)
1236      //beta_tmp = beta_uh_dry;
1237      beta_tmp = beta_uh_dry + (beta_uh - beta_uh_dry) * hfactor;
1238
1239      // Limit the gradient
1240      limit_gradient(dqv,qmin,qmax,beta_tmp);
1241
1242      for (i=0;i<3;i++)
1243        xmom_vertex_values[k3+i]=xmom_centroid_values[k]+dqv[i];
1244     
1245     
1246      //-----------------------------------
1247      // ymomentum
1248      //-----------------------------------                 
1249
1250      // Calculate the difference between vertex 0 of the auxiliary
1251      // triangle and the centroid of triangle k
1252      dq0=ymom_centroid_values[k0]-ymom_centroid_values[k];
1253     
1254      // Calculate differentials between the vertices
1255      // of the auxiliary triangle
1256      dq1=ymom_centroid_values[k1]-ymom_centroid_values[k0];
1257      dq2=ymom_centroid_values[k2]-ymom_centroid_values[k0];
1258     
1259      // Calculate the gradient of xmom on the auxiliary triangle
1260      a = dy2*dq1 - dy1*dq2;
1261      a /= area2;
1262      b = dx1*dq2 - dx2*dq1;
1263      b /= area2;
1264     
1265      // Calculate provisional jumps in stage from the centroid
1266      // of triangle k to its vertices, to be limited
1267      dqv[0]=a*dxv0+b*dyv0;
1268      dqv[1]=a*dxv1+b*dyv1;
1269      dqv[2]=a*dxv2+b*dyv2;
1270     
1271      // Now we want to find min and max of the centroid and the
1272      // vertices of the auxiliary triangle and compute jumps
1273      // from the centroid to the min and max
1274      find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
1275     
1276      //beta_tmp = beta_vh;
1277      //
1278      //if (hmin<minimum_allowed_height)
1279      //beta_tmp = beta_vh_dry;
1280      beta_tmp = beta_vh_dry + (beta_vh - beta_vh_dry) * hfactor;       
1281
1282      // Limit the gradient
1283      limit_gradient(dqv,qmin,qmax,beta_tmp);
1284     
1285      for (i=0;i<3;i++)
1286        ymom_vertex_values[k3+i]=ymom_centroid_values[k]+dqv[i];
1287       
1288    } // End number_of_boundaries <=1
1289    else{
1290
1291      //==============================================
1292      // Number of boundaries == 2
1293      //==============================================       
1294       
1295      // One internal neighbour and gradient is in direction of the neighbour's centroid
1296     
1297      // Find the only internal neighbour (k1?)
1298      for (k2=k3;k2<k3+3;k2++){
1299        // Find internal neighbour of triangle k     
1300        // k2 indexes the edges of triangle k   
1301     
1302        if (surrogate_neighbours[k2]!=k)
1303          break;
1304      }
1305     
1306      if ((k2==k3+3)) {
1307        // If we didn't find an internal neighbour
1308        PyErr_SetString(PyExc_RuntimeError, 
1309                        "shallow_water_ext.c: Internal neighbour not found");     
1310        return -1;
1311      }
1312     
1313      k1=surrogate_neighbours[k2];
1314     
1315      // The coordinates of the triangle are already (x,y).
1316      // Get centroid of the neighbour (x1,y1)
1317      coord_index=2*k1;
1318      x1=centroid_coordinates[coord_index];
1319      y1=centroid_coordinates[coord_index+1];
1320     
1321      // Compute x- and y- distances between the centroid of
1322      // triangle k and that of its neighbour
1323      dx1=x1-x; dy1=y1-y;
1324     
1325      // Set area2 as the square of the distance
1326      area2=dx1*dx1+dy1*dy1;
1327     
1328      // Set dx2=(x1-x0)/((x1-x0)^2+(y1-y0)^2)
1329      // and dy2=(y1-y0)/((x1-x0)^2+(y1-y0)^2) which
1330      // respectively correspond to the x- and y- gradients
1331      // of the conserved quantities
1332      dx2=1.0/area2;
1333      dy2=dx2*dy1;
1334      dx2*=dx1;
1335     
1336     
1337      //-----------------------------------
1338      // stage
1339      //-----------------------------------           
1340
1341      // Compute differentials
1342      dq1=stage_centroid_values[k1]-stage_centroid_values[k];
1343     
1344      // Calculate the gradient between the centroid of triangle k
1345      // and that of its neighbour
1346      a=dq1*dx2;
1347      b=dq1*dy2;
1348     
1349      // Calculate provisional vertex jumps, to be limited
1350      dqv[0]=a*dxv0+b*dyv0;
1351      dqv[1]=a*dxv1+b*dyv1;
1352      dqv[2]=a*dxv2+b*dyv2;
1353     
1354      // Now limit the jumps
1355      if (dq1>=0.0){
1356        qmin=0.0;
1357        qmax=dq1;
1358      }
1359      else{
1360        qmin=dq1;
1361        qmax=0.0;
1362      }
1363     
1364      // Limit the gradient
1365      limit_gradient(dqv,qmin,qmax,beta_w);
1366     
1367      for (i=0;i<3;i++)
1368        stage_vertex_values[k3+i]=stage_centroid_values[k]+dqv[i];
1369     
1370     
1371      //-----------------------------------
1372      // xmomentum
1373      //-----------------------------------                       
1374     
1375      // Compute differentials
1376      dq1=xmom_centroid_values[k1]-xmom_centroid_values[k];
1377     
1378      // Calculate the gradient between the centroid of triangle k
1379      // and that of its neighbour
1380      a=dq1*dx2;
1381      b=dq1*dy2;
1382     
1383      // Calculate provisional vertex jumps, to be limited
1384      dqv[0]=a*dxv0+b*dyv0;
1385      dqv[1]=a*dxv1+b*dyv1;
1386      dqv[2]=a*dxv2+b*dyv2;
1387     
1388      // Now limit the jumps
1389      if (dq1>=0.0){
1390        qmin=0.0;
1391        qmax=dq1;
1392      }
1393      else{
1394        qmin=dq1;
1395        qmax=0.0;
1396      }
1397     
1398      // Limit the gradient     
1399      limit_gradient(dqv,qmin,qmax,beta_w);
1400     
1401      for (i=0;i<3;i++)
1402        xmom_vertex_values[k3+i]=xmom_centroid_values[k]+dqv[i];
1403     
1404     
1405      //-----------------------------------
1406      // ymomentum
1407      //-----------------------------------                       
1408
1409      // Compute differentials
1410      dq1=ymom_centroid_values[k1]-ymom_centroid_values[k];
1411     
1412      // Calculate the gradient between the centroid of triangle k
1413      // and that of its neighbour
1414      a=dq1*dx2;
1415      b=dq1*dy2;
1416     
1417      // Calculate provisional vertex jumps, to be limited
1418      dqv[0]=a*dxv0+b*dyv0;
1419      dqv[1]=a*dxv1+b*dyv1;
1420      dqv[2]=a*dxv2+b*dyv2;
1421     
1422      // Now limit the jumps
1423      if (dq1>=0.0){
1424        qmin=0.0;
1425        qmax=dq1;
1426      }
1427      else{
1428        qmin=dq1;
1429        qmax=0.0;
1430      }
1431     
1432      // Limit the gradient
1433      limit_gradient(dqv,qmin,qmax,beta_w);
1434     
1435      for (i=0;i<3;i++)
1436        ymom_vertex_values[k3+i]=ymom_centroid_values[k]+dqv[i];
1437       
1438    } // else [number_of_boundaries==2]
1439  } // for k=0 to number_of_elements-1
1440 
1441  return 0;
1442}                       
1443
1444
1445PyObject *extrapolate_second_order_sw(PyObject *self, PyObject *args) {
1446  /*Compute the vertex values based on a linear reconstruction
1447    on each triangle
1448   
1449    These values are calculated as follows:
1450    1) For each triangle not adjacent to a boundary, we consider the
1451       auxiliary triangle formed by the centroids of its three
1452       neighbours.
1453    2) For each conserved quantity, we integrate around the auxiliary
1454       triangle's boundary the product of the quantity and the outward
1455       normal vector. Dividing by the triangle area gives (a,b), the
1456       average of the vector (q_x,q_y) on the auxiliary triangle.
1457       We suppose that the linear reconstruction on the original
1458       triangle has gradient (a,b).
1459    3) Provisional vertex jumps dqv[0,1,2] are computed and these are
1460       then limited by calling the functions find_qmin_and_qmax and
1461       limit_gradient
1462
1463    Python call:
1464    extrapolate_second_order_sw(domain.surrogate_neighbours,
1465                                domain.number_of_boundaries
1466                                domain.centroid_coordinates,
1467                                Stage.centroid_values
1468                                Xmom.centroid_values
1469                                Ymom.centroid_values
1470                                domain.vertex_coordinates,
1471                                Stage.vertex_values,
1472                                Xmom.vertex_values,
1473                                Ymom.vertex_values)
1474
1475    Post conditions:
1476            The vertices of each triangle have values from a
1477            limited linear reconstruction
1478            based on centroid values
1479
1480  */
1481  PyArrayObject *surrogate_neighbours,
1482    *number_of_boundaries,
1483    *centroid_coordinates,
1484    *stage_centroid_values,
1485    *xmom_centroid_values,
1486    *ymom_centroid_values,
1487    *elevation_centroid_values,
1488    *vertex_coordinates,
1489    *stage_vertex_values,
1490    *xmom_vertex_values,
1491    *ymom_vertex_values,
1492    *elevation_vertex_values;
1493 
1494  PyObject *domain;
1495
1496 
1497  double beta_w, beta_w_dry, beta_uh, beta_uh_dry, beta_vh, beta_vh_dry;   
1498  double minimum_allowed_height, epsilon;
1499  int optimise_dry_cells, number_of_elements, e;
1500 
1501  // Provisional jumps from centroids to v'tices and safety factor re limiting
1502  // by which these jumps are limited
1503  // Convert Python arguments to C
1504  if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOi",
1505                        &domain,
1506                        &surrogate_neighbours,
1507                        &number_of_boundaries,
1508                        &centroid_coordinates,
1509                        &stage_centroid_values,
1510                        &xmom_centroid_values,
1511                        &ymom_centroid_values,
1512                        &elevation_centroid_values,
1513                        &vertex_coordinates,
1514                        &stage_vertex_values,
1515                        &xmom_vertex_values,
1516                        &ymom_vertex_values,
1517                        &elevation_vertex_values,
1518                        &optimise_dry_cells)) {                 
1519                       
1520    PyErr_SetString(PyExc_RuntimeError, 
1521                    "Input arguments to extrapolate_second_order_sw failed");
1522    return NULL;
1523  }
1524
1525  // check that numpy array objects arrays are C contiguous memory
1526  CHECK_C_CONTIG(surrogate_neighbours);
1527  CHECK_C_CONTIG(number_of_boundaries);
1528  CHECK_C_CONTIG(centroid_coordinates);
1529  CHECK_C_CONTIG(stage_centroid_values);
1530  CHECK_C_CONTIG(xmom_centroid_values);
1531  CHECK_C_CONTIG(ymom_centroid_values);
1532  CHECK_C_CONTIG(elevation_centroid_values);
1533  CHECK_C_CONTIG(vertex_coordinates);
1534  CHECK_C_CONTIG(stage_vertex_values);
1535  CHECK_C_CONTIG(xmom_vertex_values);
1536  CHECK_C_CONTIG(ymom_vertex_values);
1537  CHECK_C_CONTIG(elevation_vertex_values);
1538 
1539  // Get the safety factor beta_w, set in the config.py file.
1540  // This is used in the limiting process
1541 
1542
1543  beta_w                 = get_python_double(domain,"beta_w");
1544  beta_w_dry             = get_python_double(domain,"beta_w_dry");
1545  beta_uh                = get_python_double(domain,"beta_uh");
1546  beta_uh_dry            = get_python_double(domain,"beta_uh_dry");
1547  beta_vh                = get_python_double(domain,"beta_vh");
1548  beta_vh_dry            = get_python_double(domain,"beta_vh_dry"); 
1549
1550  minimum_allowed_height = get_python_double(domain,"minimum_allowed_height");
1551  epsilon                = get_python_double(domain,"epsilon");
1552
1553  number_of_elements = stage_centroid_values -> dimensions[0]; 
1554
1555  // Call underlying computational routine
1556  e = _extrapolate_second_order_sw(number_of_elements,
1557                                   epsilon,
1558                                   minimum_allowed_height,
1559                                   beta_w,
1560                                   beta_w_dry,
1561                                   beta_uh,
1562                                   beta_uh_dry,
1563                                   beta_vh,
1564                                   beta_vh_dry,
1565                                   (long*) surrogate_neighbours -> data,
1566                                   (long*) number_of_boundaries -> data,
1567                                   (double*) centroid_coordinates -> data,
1568                                   (double*) stage_centroid_values -> data,
1569                                   (double*) xmom_centroid_values -> data,
1570                                   (double*) ymom_centroid_values -> data,
1571                                   (double*) elevation_centroid_values -> data,
1572                                   (double*) vertex_coordinates -> data,
1573                                   (double*) stage_vertex_values -> data,
1574                                   (double*) xmom_vertex_values -> data,
1575                                   (double*) ymom_vertex_values -> data,
1576                                   (double*) elevation_vertex_values -> data,
1577                                   optimise_dry_cells);
1578  if (e == -1) {
1579    // Use error string set inside computational routine
1580    return NULL;
1581  }                               
1582 
1583 
1584  return Py_BuildValue("");
1585 
1586}// extrapolate_second-order_sw
1587
1588
1589
1590
1591PyObject *rotate(PyObject *self, PyObject *args, PyObject *kwargs) {
1592  //
1593  // r = rotate(q, normal, direction=1)
1594  //
1595  // Where q is assumed to be a Float numeric array of length 3 and
1596  // normal a Float numeric array of length 2.
1597
1598  // FIXME(Ole): I don't think this is used anymore
1599
1600  PyObject *Q, *Normal;
1601  PyArrayObject *q, *r, *normal;
1602
1603  static char *argnames[] = {"q", "normal", "direction", NULL};
1604  int dimensions[1], i, direction=1;
1605  double n1, n2;
1606
1607  // Convert Python arguments to C
1608  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OO|i", argnames,
1609                                   &Q, &Normal, &direction)) {
1610    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: rotate could not parse input arguments");
1611    return NULL;
1612  } 
1613
1614  // Input checks (convert sequences into numeric arrays)
1615  q = (PyArrayObject *)
1616    PyArray_ContiguousFromObject(Q, PyArray_DOUBLE, 0, 0);
1617  normal = (PyArrayObject *)
1618    PyArray_ContiguousFromObject(Normal, PyArray_DOUBLE, 0, 0);
1619
1620
1621  if (normal -> dimensions[0] != 2) {
1622    PyErr_SetString(PyExc_RuntimeError, "Normal vector must have 2 components");
1623    return NULL;
1624  }
1625
1626  // Allocate space for return vector r (don't DECREF)
1627  dimensions[0] = 3;
1628  r = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
1629
1630  // Copy
1631  for (i=0; i<3; i++) {
1632    ((double *) (r -> data))[i] = ((double *) (q -> data))[i];
1633  }
1634
1635  // Get normal and direction
1636  n1 = ((double *) normal -> data)[0];
1637  n2 = ((double *) normal -> data)[1];
1638  if (direction == -1) n2 = -n2;
1639
1640  // Rotate
1641  _rotate((double *) r -> data, n1, n2);
1642
1643  // Release numeric arrays
1644  Py_DECREF(q);
1645  Py_DECREF(normal);
1646
1647  // Return result using PyArray to avoid memory leak
1648  return PyArray_Return(r);
1649}
1650
1651
1652// Computational function for flux computation
1653double _compute_fluxes_central(int number_of_elements,
1654                               double timestep,
1655                               double epsilon,
1656                               double H0,
1657                               double g,
1658                               long* neighbours,
1659                               long* neighbour_edges,
1660                               double* normals,
1661                               double* edgelengths, 
1662                               double* radii, 
1663                               double* areas,
1664                               long* tri_full_flag,
1665                               double* stage_edge_values,
1666                               double* xmom_edge_values,
1667                               double* ymom_edge_values,
1668                               double* bed_edge_values,
1669                               double* stage_boundary_values,
1670                               double* xmom_boundary_values,
1671                               double* ymom_boundary_values,
1672                               double* stage_explicit_update,
1673                               double* xmom_explicit_update,
1674                               double* ymom_explicit_update,
1675                               long* already_computed_flux,
1676                               double* max_speed_array,
1677                               int optimise_dry_cells) {
1678                               
1679  // Local variables
1680  double max_speed, length, area, zl, zr;
1681  int k, i, m, n;
1682  int ki, nm=0, ki2; // Index shorthands
1683 
1684  // Workspace (making them static actually made function slightly slower (Ole)) 
1685  double ql[3], qr[3], edgeflux[3]; // Work array for summing up fluxes
1686
1687  static long call=1; // Static local variable flagging already computed flux
1688                               
1689       
1690  // Start computation
1691  call++; // Flag 'id' of flux calculation for this timestep
1692 
1693  // Set explicit_update to zero for all conserved_quantities.
1694  // This assumes compute_fluxes called before forcing terms
1695  for (k=0; k<number_of_elements; k++) {
1696    stage_explicit_update[k]=0.0;
1697    xmom_explicit_update[k]=0.0;
1698    ymom_explicit_update[k]=0.0; 
1699  }
1700
1701  // For all triangles
1702  for (k=0; k<number_of_elements; k++) {
1703   
1704    // Loop through neighbours and compute edge flux for each 
1705    for (i=0; i<3; i++) {
1706      ki = k*3+i; // Linear index (triangle k, edge i)
1707     
1708      if (already_computed_flux[ki] == call)
1709        // We've already computed the flux across this edge
1710        continue;
1711       
1712       
1713      ql[0] = stage_edge_values[ki];
1714      ql[1] = xmom_edge_values[ki];
1715      ql[2] = ymom_edge_values[ki];
1716      zl = bed_edge_values[ki];
1717
1718      // Quantities at neighbour on nearest face
1719      n = neighbours[ki];
1720      if (n < 0) {
1721        // Neighbour is a boundary condition
1722        m = -n-1; // Convert negative flag to boundary index
1723       
1724        qr[0] = stage_boundary_values[m];
1725        qr[1] = xmom_boundary_values[m];
1726        qr[2] = ymom_boundary_values[m];
1727        zr = zl; // Extend bed elevation to boundary
1728      } else {
1729        // Neighbour is a real element
1730        m = neighbour_edges[ki];
1731        nm = n*3+m; // Linear index (triangle n, edge m)
1732       
1733        qr[0] = stage_edge_values[nm];
1734        qr[1] = xmom_edge_values[nm];
1735        qr[2] = ymom_edge_values[nm];
1736        zr = bed_edge_values[nm];
1737      }
1738     
1739     
1740      if (optimise_dry_cells) {     
1741        // Check if flux calculation is necessary across this edge
1742        // This check will exclude dry cells.
1743        // This will also optimise cases where zl != zr as
1744        // long as both are dry
1745
1746        if ( fabs(ql[0] - zl) < epsilon && 
1747             fabs(qr[0] - zr) < epsilon ) {
1748          // Cell boundary is dry
1749         
1750          already_computed_flux[ki] = call; // #k Done 
1751          if (n>=0)
1752            already_computed_flux[nm] = call; // #n Done
1753       
1754          max_speed = 0.0;
1755          continue;
1756        }
1757      }
1758   
1759     
1760      // Outward pointing normal vector (domain.normals[k, 2*i:2*i+2])
1761      ki2 = 2*ki; //k*6 + i*2
1762
1763      // Edge flux computation (triangle k, edge i)
1764      _flux_function_central(ql, qr, zl, zr,
1765                             normals[ki2], normals[ki2+1],
1766                             epsilon, H0, g,
1767                             edgeflux, &max_speed);
1768     
1769     
1770      // Multiply edgeflux by edgelength
1771      length = edgelengths[ki];
1772      edgeflux[0] *= length;           
1773      edgeflux[1] *= length;           
1774      edgeflux[2] *= length;                       
1775     
1776     
1777      // Update triangle k with flux from edge i
1778      stage_explicit_update[k] -= edgeflux[0];
1779      xmom_explicit_update[k] -= edgeflux[1];
1780      ymom_explicit_update[k] -= edgeflux[2];
1781     
1782      already_computed_flux[ki] = call; // #k Done
1783     
1784     
1785      // Update neighbour n with same flux but reversed sign
1786      if (n>=0) {
1787        stage_explicit_update[n] += edgeflux[0];
1788        xmom_explicit_update[n] += edgeflux[1];
1789        ymom_explicit_update[n] += edgeflux[2];
1790       
1791        already_computed_flux[nm] = call; // #n Done
1792      }
1793
1794
1795      // Update timestep based on edge i and possibly neighbour n
1796      if (tri_full_flag[k] == 1) {
1797        if (max_speed > epsilon) {
1798       
1799          // Apply CFL condition for triangles joining this edge (triangle k and triangle n)
1800         
1801          // CFL for triangle k
1802          timestep = min(timestep, radii[k]/max_speed);
1803         
1804          if (n>=0) 
1805            // Apply CFL condition for neigbour n (which is on the ith edge of triangle k)
1806            timestep = min(timestep, radii[n]/max_speed);
1807         
1808          // Ted Rigby's suggested less conservative version
1809          //if (n>=0) {             
1810          //  timestep = min(timestep, (radii[k]+radii[n])/max_speed);
1811          //} else {
1812          //  timestep = min(timestep, radii[k]/max_speed);
1813          // }
1814        }
1815      }
1816     
1817    } // End edge i (and neighbour n)
1818   
1819   
1820    // Normalise triangle k by area and store for when all conserved
1821    // quantities get updated
1822    area = areas[k];
1823    stage_explicit_update[k] /= area;
1824    xmom_explicit_update[k] /= area;
1825    ymom_explicit_update[k] /= area;
1826   
1827   
1828    // Keep track of maximal speeds
1829    max_speed_array[k] = max_speed;
1830   
1831  } // End triangle k
1832       
1833                               
1834                               
1835  return timestep;
1836}
1837
1838//=========================================================================
1839// Python Glue
1840//=========================================================================
1841
1842PyObject *compute_fluxes_ext_central_new(PyObject *self, PyObject *args) {
1843  /*Compute all fluxes and the timestep suitable for all volumes
1844    in domain.
1845
1846    Compute total flux for each conserved quantity using "flux_function_central"
1847
1848    Fluxes across each edge are scaled by edgelengths and summed up
1849    Resulting flux is then scaled by area and stored in
1850    explicit_update for each of the three conserved quantities
1851    stage, xmomentum and ymomentum
1852
1853    The maximal allowable speed computed by the flux_function for each volume
1854    is converted to a timestep that must not be exceeded. The minimum of
1855    those is computed as the next overall timestep.
1856
1857    Python call:
1858    timestep = compute_fluxes(timestep, domain, stage, xmom, ymom, bed)
1859
1860
1861    Post conditions:
1862      domain.explicit_update is reset to computed flux values
1863      returns timestep which is the largest step satisfying all volumes.
1864
1865
1866  */
1867
1868    PyObject
1869        *domain,
1870        *stage, 
1871        *xmom, 
1872        *ymom, 
1873        *bed;
1874
1875    PyArrayObject
1876        *neighbours, 
1877        *neighbour_edges,
1878        *normals, 
1879        *edgelengths, 
1880        *radii, 
1881        *areas,
1882        *tri_full_flag,
1883        *stage_edge_values,
1884        *xmom_edge_values,
1885        *ymom_edge_values,
1886        *bed_edge_values,
1887        *stage_boundary_values,
1888        *xmom_boundary_values,
1889        *ymom_boundary_values,
1890        *stage_explicit_update,
1891        *xmom_explicit_update,
1892        *ymom_explicit_update,
1893        *already_computed_flux, //Tracks whether the flux across an edge has already been computed
1894        *max_speed_array; //Keeps track of max speeds for each triangle
1895
1896   
1897    double timestep, epsilon, H0, g;
1898    int optimise_dry_cells;
1899   
1900    // Convert Python arguments to C
1901    if (!PyArg_ParseTuple(args, "dOOOO", &timestep, &domain, &stage, &xmom, &ymom, &bed )) {
1902        PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
1903        return NULL;
1904    }
1905
1906    epsilon           = get_python_double(domain,"epsilon");
1907    H0                = get_python_double(domain,"H0");
1908    g                 = get_python_double(domain,"g");
1909    optimise_dry_cells = get_python_integer(domain,"optimse_dry_cells");
1910   
1911    neighbours        = get_consecutive_array(domain, "neighbours");
1912    neighbour_edges   = get_consecutive_array(domain, "neighbour_edges"); 
1913    normals           = get_consecutive_array(domain, "normals");
1914    edgelengths       = get_consecutive_array(domain, "edge_lengths");   
1915    radii             = get_consecutive_array(domain, "radii");   
1916    areas             = get_consecutive_array(domain, "areas");   
1917    tri_full_flag     = get_consecutive_array(domain, "tri_full_flag");
1918    already_computed_flux  = get_consecutive_array(domain, "already_computed_flux");
1919    max_speed_array   = get_consecutive_array(domain, "max_speed");
1920   
1921    stage_edge_values = get_consecutive_array(stage, "edge_values");   
1922    xmom_edge_values  = get_consecutive_array(xmom, "edge_values");   
1923    ymom_edge_values  = get_consecutive_array(ymom, "edge_values"); 
1924    bed_edge_values   = get_consecutive_array(bed, "edge_values");   
1925
1926    stage_boundary_values = get_consecutive_array(stage, "boundary_values");   
1927    xmom_boundary_values  = get_consecutive_array(xmom, "boundary_values");   
1928    ymom_boundary_values  = get_consecutive_array(ymom, "boundary_values"); 
1929
1930    stage_explicit_update = get_consecutive_array(stage, "explicit_update");   
1931    xmom_explicit_update  = get_consecutive_array(xmom, "explicit_update");   
1932    ymom_explicit_update  = get_consecutive_array(ymom, "explicit_update"); 
1933
1934
1935  int number_of_elements = stage_edge_values -> dimensions[0];
1936
1937  // Call underlying flux computation routine and update
1938  // the explicit update arrays
1939  timestep = _compute_fluxes_central(number_of_elements,
1940                                     timestep,
1941                                     epsilon,
1942                                     H0,
1943                                     g,
1944                                     (long*) neighbours -> data,
1945                                     (long*) neighbour_edges -> data,
1946                                     (double*) normals -> data,
1947                                     (double*) edgelengths -> data, 
1948                                     (double*) radii -> data, 
1949                                     (double*) areas -> data,
1950                                     (long*) tri_full_flag -> data,
1951                                     (double*) stage_edge_values -> data,
1952                                     (double*) xmom_edge_values -> data,
1953                                     (double*) ymom_edge_values -> data,
1954                                     (double*) bed_edge_values -> data,
1955                                     (double*) stage_boundary_values -> data,
1956                                     (double*) xmom_boundary_values -> data,
1957                                     (double*) ymom_boundary_values -> data,
1958                                     (double*) stage_explicit_update -> data,
1959                                     (double*) xmom_explicit_update -> data,
1960                                     (double*) ymom_explicit_update -> data,
1961                                     (long*) already_computed_flux -> data,
1962                                     (double*) max_speed_array -> data,
1963                                     optimise_dry_cells);
1964
1965  Py_DECREF(neighbours);
1966  Py_DECREF(neighbour_edges);
1967  Py_DECREF(normals);
1968  Py_DECREF(edgelengths);
1969  Py_DECREF(radii);
1970  Py_DECREF(areas);
1971  Py_DECREF(tri_full_flag);
1972  Py_DECREF(already_computed_flux);
1973  Py_DECREF(max_speed_array);
1974  Py_DECREF(stage_edge_values);
1975  Py_DECREF(xmom_edge_values);
1976  Py_DECREF(ymom_edge_values);
1977  Py_DECREF(bed_edge_values);
1978  Py_DECREF(stage_boundary_values);
1979  Py_DECREF(xmom_boundary_values);
1980  Py_DECREF(ymom_boundary_values);
1981  Py_DECREF(stage_explicit_update);
1982  Py_DECREF(xmom_explicit_update);
1983  Py_DECREF(ymom_explicit_update);
1984
1985 
1986  // Return updated flux timestep
1987  return Py_BuildValue("d", timestep);
1988}
1989
1990
1991
1992
1993
1994
1995PyObject *compute_fluxes_ext_central(PyObject *self, PyObject *args) {
1996  /*Compute all fluxes and the timestep suitable for all volumes
1997    in domain.
1998
1999    Compute total flux for each conserved quantity using "flux_function_central"
2000
2001    Fluxes across each edge are scaled by edgelengths and summed up
2002    Resulting flux is then scaled by area and stored in
2003    explicit_update for each of the three conserved quantities
2004    stage, xmomentum and ymomentum
2005
2006    The maximal allowable speed computed by the flux_function for each volume
2007    is converted to a timestep that must not be exceeded. The minimum of
2008    those is computed as the next overall timestep.
2009
2010    Python call:
2011    domain.timestep = compute_fluxes(timestep,
2012                                     domain.epsilon,
2013                                     domain.H0,
2014                                     domain.g,
2015                                     domain.neighbours,
2016                                     domain.neighbour_edges,
2017                                     domain.normals,
2018                                     domain.edgelengths,
2019                                     domain.radii,
2020                                     domain.areas,
2021                                     tri_full_flag,
2022                                     Stage.edge_values,
2023                                     Xmom.edge_values,
2024                                     Ymom.edge_values,
2025                                     Bed.edge_values,
2026                                     Stage.boundary_values,
2027                                     Xmom.boundary_values,
2028                                     Ymom.boundary_values,
2029                                     Stage.explicit_update,
2030                                     Xmom.explicit_update,
2031                                     Ymom.explicit_update,
2032                                     already_computed_flux,
2033                                     optimise_dry_cells)                                     
2034
2035
2036    Post conditions:
2037      domain.explicit_update is reset to computed flux values
2038      domain.timestep is set to the largest step satisfying all volumes.
2039
2040
2041  */
2042
2043
2044  PyArrayObject *neighbours, *neighbour_edges,
2045    *normals, *edgelengths, *radii, *areas,
2046    *tri_full_flag,
2047    *stage_edge_values,
2048    *xmom_edge_values,
2049    *ymom_edge_values,
2050    *bed_edge_values,
2051    *stage_boundary_values,
2052    *xmom_boundary_values,
2053    *ymom_boundary_values,
2054    *stage_explicit_update,
2055    *xmom_explicit_update,
2056    *ymom_explicit_update,
2057    *already_computed_flux, //Tracks whether the flux across an edge has already been computed
2058    *max_speed_array; //Keeps track of max speeds for each triangle
2059
2060   
2061  double timestep, epsilon, H0, g;
2062  int optimise_dry_cells;
2063   
2064  // Convert Python arguments to C
2065  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOi",
2066                        &timestep,
2067                        &epsilon,
2068                        &H0,
2069                        &g,
2070                        &neighbours,
2071                        &neighbour_edges,
2072                        &normals,
2073                        &edgelengths, &radii, &areas,
2074                        &tri_full_flag,
2075                        &stage_edge_values,
2076                        &xmom_edge_values,
2077                        &ymom_edge_values,
2078                        &bed_edge_values,
2079                        &stage_boundary_values,
2080                        &xmom_boundary_values,
2081                        &ymom_boundary_values,
2082                        &stage_explicit_update,
2083                        &xmom_explicit_update,
2084                        &ymom_explicit_update,
2085                        &already_computed_flux,
2086                        &max_speed_array,
2087                        &optimise_dry_cells)) {
2088    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
2089    return NULL;
2090  }
2091
2092  // check that numpy array objects arrays are C contiguous memory
2093  CHECK_C_CONTIG(neighbours);
2094  CHECK_C_CONTIG(neighbour_edges);
2095  CHECK_C_CONTIG(normals);
2096  CHECK_C_CONTIG(edgelengths);
2097  CHECK_C_CONTIG(radii);
2098  CHECK_C_CONTIG(areas);
2099  CHECK_C_CONTIG(tri_full_flag);
2100  CHECK_C_CONTIG(stage_edge_values);
2101  CHECK_C_CONTIG(xmom_edge_values);
2102  CHECK_C_CONTIG(ymom_edge_values);
2103  CHECK_C_CONTIG(bed_edge_values);
2104  CHECK_C_CONTIG(stage_boundary_values);
2105  CHECK_C_CONTIG(xmom_boundary_values);
2106  CHECK_C_CONTIG(ymom_boundary_values);
2107  CHECK_C_CONTIG(stage_explicit_update);
2108  CHECK_C_CONTIG(xmom_explicit_update);
2109  CHECK_C_CONTIG(ymom_explicit_update);
2110  CHECK_C_CONTIG(already_computed_flux);
2111  CHECK_C_CONTIG(max_speed_array);
2112
2113  int number_of_elements = stage_edge_values -> dimensions[0];
2114
2115  // Call underlying flux computation routine and update
2116  // the explicit update arrays
2117  timestep = _compute_fluxes_central(number_of_elements,
2118                                     timestep,
2119                                     epsilon,
2120                                     H0,
2121                                     g,
2122                                     (long*) neighbours -> data,
2123                                     (long*) neighbour_edges -> data,
2124                                     (double*) normals -> data,
2125                                     (double*) edgelengths -> data, 
2126                                     (double*) radii -> data, 
2127                                     (double*) areas -> data,
2128                                     (long*) tri_full_flag -> data,
2129                                     (double*) stage_edge_values -> data,
2130                                     (double*) xmom_edge_values -> data,
2131                                     (double*) ymom_edge_values -> data,
2132                                     (double*) bed_edge_values -> data,
2133                                     (double*) stage_boundary_values -> data,
2134                                     (double*) xmom_boundary_values -> data,
2135                                     (double*) ymom_boundary_values -> data,
2136                                     (double*) stage_explicit_update -> data,
2137                                     (double*) xmom_explicit_update -> data,
2138                                     (double*) ymom_explicit_update -> data,
2139                                     (long*) already_computed_flux -> data,
2140                                     (double*) max_speed_array -> data,
2141                                     optimise_dry_cells);
2142 
2143  // Return updated flux timestep
2144  return Py_BuildValue("d", timestep);
2145}
2146
2147
2148
2149
2150
2151PyObject *compute_fluxes_ext_kinetic(PyObject *self, PyObject *args) {
2152  /*Compute all fluxes and the timestep suitable for all volumes
2153    in domain.
2154
2155    THIS IS AN EXPERIMENTAL FUNCTION - NORMALLY flux_function_central IS USED.
2156   
2157    Compute total flux for each conserved quantity using "flux_function_kinetic"
2158
2159    Fluxes across each edge are scaled by edgelengths and summed up
2160    Resulting flux is then scaled by area and stored in
2161    explicit_update for each of the three conserved quantities
2162    stage, xmomentum and ymomentum
2163
2164    The maximal allowable speed computed by the flux_function for each volume
2165    is converted to a timestep that must not be exceeded. The minimum of
2166    those is computed as the next overall timestep.
2167
2168    Python call:
2169    domain.timestep = compute_fluxes(timestep,
2170                                     domain.epsilon,
2171                                     domain.H0,
2172                                     domain.g,
2173                                     domain.neighbours,
2174                                     domain.neighbour_edges,
2175                                     domain.normals,
2176                                     domain.edgelengths,
2177                                     domain.radii,
2178                                     domain.areas,
2179                                     Stage.edge_values,
2180                                     Xmom.edge_values,
2181                                     Ymom.edge_values,
2182                                     Bed.edge_values,
2183                                     Stage.boundary_values,
2184                                     Xmom.boundary_values,
2185                                     Ymom.boundary_values,
2186                                     Stage.explicit_update,
2187                                     Xmom.explicit_update,
2188                                     Ymom.explicit_update,
2189                                     already_computed_flux)
2190
2191
2192    Post conditions:
2193      domain.explicit_update is reset to computed flux values
2194      domain.timestep is set to the largest step satisfying all volumes.
2195
2196
2197  */
2198
2199
2200  PyArrayObject *neighbours, *neighbour_edges,
2201    *normals, *edgelengths, *radii, *areas,
2202    *tri_full_flag,
2203    *stage_edge_values,
2204    *xmom_edge_values,
2205    *ymom_edge_values,
2206    *bed_edge_values,
2207    *stage_boundary_values,
2208    *xmom_boundary_values,
2209    *ymom_boundary_values,
2210    *stage_explicit_update,
2211    *xmom_explicit_update,
2212    *ymom_explicit_update,
2213    *already_computed_flux; // Tracks whether the flux across an edge has already been computed
2214
2215
2216  // Local variables
2217  double timestep, max_speed, epsilon, g, H0;
2218  double normal[2], ql[3], qr[3], zl, zr;
2219  double edgeflux[3]; //Work arrays for summing up fluxes
2220
2221  int number_of_elements, k, i, m, n;
2222  int ki, nm=0, ki2; //Index shorthands
2223  static long call=1;
2224
2225
2226  // Convert Python arguments to C
2227  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOO",
2228                        &timestep,
2229                        &epsilon,
2230                        &H0,
2231                        &g,
2232                        &neighbours,
2233                        &neighbour_edges,
2234                        &normals,
2235                        &edgelengths, &radii, &areas,
2236                        &tri_full_flag,
2237                        &stage_edge_values,
2238                        &xmom_edge_values,
2239                        &ymom_edge_values,
2240                        &bed_edge_values,
2241                        &stage_boundary_values,
2242                        &xmom_boundary_values,
2243                        &ymom_boundary_values,
2244                        &stage_explicit_update,
2245                        &xmom_explicit_update,
2246                        &ymom_explicit_update,
2247                        &already_computed_flux)) {
2248    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
2249    return NULL;
2250  }
2251  number_of_elements = stage_edge_values -> dimensions[0];
2252  call++;//a static local variable to which already_computed_flux is compared
2253  //set explicit_update to zero for all conserved_quantities.
2254  //This assumes compute_fluxes called before forcing terms
2255  for (k=0; k<number_of_elements; k++) {
2256    ((double *) stage_explicit_update -> data)[k]=0.0;
2257    ((double *) xmom_explicit_update -> data)[k]=0.0;
2258    ((double *) ymom_explicit_update -> data)[k]=0.0;
2259  }
2260  //Loop through neighbours and compute edge flux for each
2261  for (k=0; k<number_of_elements; k++) {
2262    for (i=0; i<3; i++) {
2263      ki = k*3+i;
2264      if (((long *) already_computed_flux->data)[ki]==call)//we've already computed the flux across this edge
2265        continue;
2266      ql[0] = ((double *) stage_edge_values -> data)[ki];
2267      ql[1] = ((double *) xmom_edge_values -> data)[ki];
2268      ql[2] = ((double *) ymom_edge_values -> data)[ki];
2269      zl =    ((double *) bed_edge_values -> data)[ki];
2270
2271      //Quantities at neighbour on nearest face
2272      n = ((long *) neighbours -> data)[ki];
2273      if (n < 0) {
2274        m = -n-1; //Convert negative flag to index
2275        qr[0] = ((double *) stage_boundary_values -> data)[m];
2276        qr[1] = ((double *) xmom_boundary_values -> data)[m];
2277        qr[2] = ((double *) ymom_boundary_values -> data)[m];
2278        zr = zl; //Extend bed elevation to boundary
2279      } else {
2280        m = ((long *) neighbour_edges -> data)[ki];
2281        nm = n*3+m;
2282        qr[0] = ((double *) stage_edge_values -> data)[nm];
2283        qr[1] = ((double *) xmom_edge_values -> data)[nm];
2284        qr[2] = ((double *) ymom_edge_values -> data)[nm];
2285        zr =    ((double *) bed_edge_values -> data)[nm];
2286      }
2287      // Outward pointing normal vector
2288      // normal = domain.normals[k, 2*i:2*i+2]
2289      ki2 = 2*ki; //k*6 + i*2
2290      normal[0] = ((double *) normals -> data)[ki2];
2291      normal[1] = ((double *) normals -> data)[ki2+1];
2292      //Edge flux computation
2293      flux_function_kinetic(ql, qr, zl, zr,
2294                    normal[0], normal[1],
2295                    epsilon, H0, g,
2296                    edgeflux, &max_speed);
2297      //update triangle k
2298      ((long *) already_computed_flux->data)[ki]=call;
2299      ((double *) stage_explicit_update -> data)[k] -= edgeflux[0]*((double *) edgelengths -> data)[ki];
2300      ((double *) xmom_explicit_update -> data)[k] -= edgeflux[1]*((double *) edgelengths -> data)[ki];
2301      ((double *) ymom_explicit_update -> data)[k] -= edgeflux[2]*((double *) edgelengths -> data)[ki];
2302      //update the neighbour n
2303      if (n>=0){
2304        ((long *) already_computed_flux->data)[nm]=call;
2305        ((double *) stage_explicit_update -> data)[n] += edgeflux[0]*((double *) edgelengths -> data)[nm];
2306        ((double *) xmom_explicit_update -> data)[n] += edgeflux[1]*((double *) edgelengths -> data)[nm];
2307        ((double *) ymom_explicit_update -> data)[n] += edgeflux[2]*((double *) edgelengths -> data)[nm];
2308      }
2309      ///for (j=0; j<3; j++) {
2310        ///flux[j] -= edgeflux[j]*((double *) edgelengths -> data)[ki];
2311        ///}
2312        //Update timestep
2313        //timestep = min(timestep, domain.radii[k]/max_speed)
2314        //FIXME: SR Add parameter for CFL condition
2315    if ( ((long *) tri_full_flag -> data)[k] == 1) {
2316            if (max_speed > epsilon) {
2317                timestep = min(timestep, ((double *) radii -> data)[k]/max_speed);
2318                //maxspeed in flux_function is calculated as max(|u+a|,|u-a|)
2319                if (n>=0)
2320                    timestep = min(timestep, ((double *) radii -> data)[n]/max_speed);
2321            }
2322    }
2323    } // end for i
2324    //Normalise by area and store for when all conserved
2325    //quantities get updated
2326    ((double *) stage_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
2327    ((double *) xmom_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
2328    ((double *) ymom_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
2329  } //end for k
2330  return Py_BuildValue("d", timestep);
2331}
2332
2333PyObject *protect(PyObject *self, PyObject *args) {
2334  //
2335  //    protect(minimum_allowed_height, maximum_allowed_speed, wc, zc, xmomc, ymomc)
2336
2337
2338  PyArrayObject
2339  *wc,            //Stage at centroids
2340  *zc,            //Elevation at centroids
2341  *xmomc,         //Momentums at centroids
2342  *ymomc;
2343
2344
2345  int N;
2346  double minimum_allowed_height, maximum_allowed_speed, epsilon;
2347
2348  // Convert Python arguments to C
2349  if (!PyArg_ParseTuple(args, "dddOOOO",
2350                        &minimum_allowed_height,
2351                        &maximum_allowed_speed,
2352                        &epsilon,
2353                        &wc, &zc, &xmomc, &ymomc)) {
2354    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: protect could not parse input arguments");
2355    return NULL;
2356  } 
2357
2358  N = wc -> dimensions[0];
2359
2360  _protect(N,
2361           minimum_allowed_height,
2362           maximum_allowed_speed,
2363           epsilon,
2364           (double*) wc -> data,
2365           (double*) zc -> data,
2366           (double*) xmomc -> data,
2367           (double*) ymomc -> data);
2368
2369  return Py_BuildValue("");
2370}
2371
2372
2373
2374PyObject *balance_deep_and_shallow(PyObject *self, PyObject *args) {
2375  // Compute linear combination between stage as computed by
2376  // gradient-limiters limiting using w, and stage computed by
2377  // gradient-limiters limiting using h (h-limiter).
2378  // The former takes precedence when heights are large compared to the
2379  // bed slope while the latter takes precedence when heights are
2380  // relatively small.  Anything in between is computed as a balanced
2381  // linear combination in order to avoid numerical disturbances which
2382  // would otherwise appear as a result of hard switching between
2383  // modes.
2384  //
2385  //    balance_deep_and_shallow(wc, zc, wv, zv,
2386  //                             xmomc, ymomc, xmomv, ymomv)
2387
2388
2389  PyArrayObject
2390    *wc,            //Stage at centroids
2391    *zc,            //Elevation at centroids
2392    *wv,            //Stage at vertices
2393    *zv,            //Elevation at vertices
2394    *hvbar,         //h-Limited depths at vertices
2395    *xmomc,         //Momentums at centroids and vertices
2396    *ymomc,
2397    *xmomv,
2398    *ymomv;
2399 
2400  PyObject *domain, *Tmp;
2401   
2402  double alpha_balance = 2.0;
2403  double H0;
2404
2405  int N, tight_slope_limiters, use_centroid_velocities; //, err;
2406
2407  // Convert Python arguments to C
2408  if (!PyArg_ParseTuple(args, "OOOOOOOOOO",
2409                        &domain,
2410                        &wc, &zc, 
2411                        &wv, &zv, &hvbar,
2412                        &xmomc, &ymomc, &xmomv, &ymomv)) {
2413    PyErr_SetString(PyExc_RuntimeError, 
2414                    "shallow_water_ext.c: balance_deep_and_shallow could not parse input arguments");
2415    return NULL;
2416  } 
2417         
2418         
2419  // FIXME (Ole): I tested this without GetAttrString and got time down
2420  // marginally from 4.0s to 3.8s. Consider passing everything in
2421  // through ParseTuple and profile.
2422 
2423  // Pull out parameters
2424  Tmp = PyObject_GetAttrString(domain, "alpha_balance");
2425  if (!Tmp) {
2426    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not obtain object alpha_balance from domain");
2427    return NULL;
2428  } 
2429  alpha_balance = PyFloat_AsDouble(Tmp);
2430  Py_DECREF(Tmp);
2431
2432 
2433  Tmp = PyObject_GetAttrString(domain, "H0");
2434  if (!Tmp) {
2435    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not obtain object H0 from domain");
2436    return NULL;
2437  } 
2438  H0 = PyFloat_AsDouble(Tmp);
2439  Py_DECREF(Tmp);
2440
2441 
2442  Tmp = PyObject_GetAttrString(domain, "tight_slope_limiters");
2443  if (!Tmp) {
2444    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not obtain object tight_slope_limiters from domain");
2445    return NULL;
2446  } 
2447  tight_slope_limiters = PyInt_AsLong(Tmp);
2448  Py_DECREF(Tmp);
2449 
2450 
2451  Tmp = PyObject_GetAttrString(domain, "use_centroid_velocities");
2452  if (!Tmp) {
2453    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not obtain object use_centroid_velocities from domain");
2454    return NULL;
2455  } 
2456  use_centroid_velocities = PyInt_AsLong(Tmp);
2457  Py_DECREF(Tmp);
2458 
2459 
2460 
2461  N = wc -> dimensions[0];
2462  _balance_deep_and_shallow(N,
2463                            (double*) wc -> data,
2464                            (double*) zc -> data,
2465                            (double*) wv -> data,
2466                            (double*) zv -> data,
2467                            (double*) hvbar -> data,
2468                            (double*) xmomc -> data,
2469                            (double*) ymomc -> data,
2470                            (double*) xmomv -> data,
2471                            (double*) ymomv -> data,
2472                            H0,
2473                            (int) tight_slope_limiters,
2474                            (int) use_centroid_velocities,                         
2475                            alpha_balance);
2476
2477
2478  return Py_BuildValue("");
2479}
2480
2481
2482
2483
2484PyObject *assign_windfield_values(PyObject *self, PyObject *args) {
2485  //
2486  //      assign_windfield_values(xmom_update, ymom_update,
2487  //                              s_vec, phi_vec, self.const)
2488
2489
2490
2491  PyArrayObject   //(one element per triangle)
2492  *s_vec,         //Speeds
2493  *phi_vec,       //Bearings
2494  *xmom_update,   //Momentum updates
2495  *ymom_update;
2496
2497
2498  int N;
2499  double cw;
2500
2501  // Convert Python arguments to C
2502  if (!PyArg_ParseTuple(args, "OOOOd",
2503                        &xmom_update,
2504                        &ymom_update,
2505                        &s_vec, &phi_vec,
2506                        &cw)) {
2507    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: assign_windfield_values could not parse input arguments");
2508    return NULL;
2509  } 
2510                       
2511
2512  N = xmom_update -> dimensions[0];
2513
2514  _assign_wind_field_values(N,
2515           (double*) xmom_update -> data,
2516           (double*) ymom_update -> data,
2517           (double*) s_vec -> data,
2518           (double*) phi_vec -> data,
2519           cw);
2520
2521  return Py_BuildValue("");
2522}
2523
2524
2525
2526
2527//-------------------------------
2528// Method table for python module
2529//-------------------------------
2530static struct PyMethodDef MethodTable[] = {
2531  /* The cast of the function is necessary since PyCFunction values
2532   * only take two PyObject* parameters, and rotate() takes
2533   * three.
2534   */
2535
2536  {"rotate", (PyCFunction)rotate, METH_VARARGS | METH_KEYWORDS, "Print out"},
2537  {"extrapolate_second_order_sw", extrapolate_second_order_sw, METH_VARARGS, "Print out"},
2538  {"compute_fluxes_ext_central", compute_fluxes_ext_central, METH_VARARGS, "Print out"},
2539  {"compute_fluxes_ext_central_new", compute_fluxes_ext_central_new, METH_VARARGS, "Print out"},
2540  {"compute_fluxes_ext_kinetic", compute_fluxes_ext_kinetic, METH_VARARGS, "Print out"},
2541  {"gravity", gravity, METH_VARARGS, "Print out"},
2542  {"manning_friction", manning_friction, METH_VARARGS, "Print out"},
2543  {"flux_function_central", flux_function_central, METH_VARARGS, "Print out"}, 
2544  {"balance_deep_and_shallow", balance_deep_and_shallow,
2545   METH_VARARGS, "Print out"},
2546  {"protect", protect, METH_VARARGS | METH_KEYWORDS, "Print out"},
2547  {"assign_windfield_values", assign_windfield_values,
2548   METH_VARARGS | METH_KEYWORDS, "Print out"},
2549  {NULL, NULL}
2550};
2551
2552// Module initialisation
2553void initshallow_water_ext(void){
2554  Py_InitModule("shallow_water_ext", MethodTable);
2555
2556  import_array(); // Necessary for handling of NumPY structures
2557}
Note: See TracBrowser for help on using the repository browser.