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

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

Refactored options, so that the use of centroid_velocities is a separate
option to tight_slope_limiters.

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