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

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

Rolled back changes from changeset:5297 regarding ticket:262.
All tests and validations passes.

File size: 76.6 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; 
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 uc, vc; // Centroid speeds
561
562  // Compute linear combination between w-limited stages and
563  // h-limited stages close to the bed elevation.
564
565  for (k=0; k<N; k++) {
566    // Compute maximal variation in bed elevation
567    // This quantitiy is
568    //     dz = max_i abs(z_i - z_c)
569    // and it is independent of dimension
570    // In the 1d case zc = (z0+z1)/2
571    // In the 2d case zc = (z0+z1+z2)/3
572
573    k3 = 3*k;
574    hc_k = wc[k] - zc[k]; // Centroid value at triangle k
575   
576    dz = 0.0;
577    if (tight_slope_limiters == 0) {     
578      // FIXME: Try with this one precomputed
579      for (i=0; i<3; i++) {
580        dz = max(dz, fabs(zv[k3+i]-zc[k]));
581      }
582    }
583
584    // Calculate depth at vertices (possibly negative here!)
585    hv[0] = wv[k3] -   zv[k3];
586    hv[1] = wv[k3+1] - zv[k3+1];
587    hv[2] = wv[k3+2] - zv[k3+2];       
588   
589    // Calculate minimal depth across all three vertices
590    hmin = min(hv[0], min(hv[1], hv[2]));
591
592    //if (hmin < 0.0 ) {
593    //  printf("hmin = %f\n",hmin);
594    //}
595   
596
597    // Create alpha in [0,1], where alpha==0 means using the h-limited
598    // stage and alpha==1 means using the w-limited stage as
599    // computed by the gradient limiter (both 1st or 2nd order)
600    if (tight_slope_limiters == 0) {     
601      // If hmin > dz/alpha_balance then alpha = 1 and the bed will have no
602      // effect
603      // If hmin < 0 then alpha = 0 reverting to constant height above bed.
604      // The parameter alpha_balance==2 by default
605
606     
607      if (dz > 0.0) {
608        alpha = max( min( alpha_balance*hmin/dz, 1.0), 0.0 );     
609      } else {
610        alpha = 1.0;  // Flat bed
611      }
612      //printf("Using old style limiter\n");
613     
614    } else {
615
616      // Tight Slope Limiter (2007)
617   
618      // Make alpha as large as possible but still ensure that
619      // final depth is positive and that velocities at vertices
620      // are controlled
621   
622      if (hmin < H0) {
623        alpha = 1.0;
624        for (i=0; i<3; i++) {
625       
626          // FIXME (Ole): Simplify (remove) when (if) hvbar gets retired
627          if (beta_h > epsilon) {
628            hc_k = hvbar[k3+i]; // Depth to be used at vertices
629          }
630         
631          h_diff = hc_k - hv[i];         
632          if (h_diff <= 0) {
633            // Deep water triangle is further away from bed than
634            // shallow water (hbar < h). Any alpha will do
635         
636          } else { 
637            // Denominator is positive which means that we need some of the
638            // h-limited stage.
639           
640            alpha = min(alpha, (hc_k - H0)/h_diff);         
641          }
642        }
643
644        // Ensure alpha in [0,1]
645        if (alpha>1.0) alpha=1.0;
646        if (alpha<0.0) alpha=0.0;
647       
648      } else {
649        // Use w-limited stage exclusively in deeper water.
650        alpha = 1.0;       
651      }
652    }
653
654
655    //  Let
656    //
657    //    wvi be the w-limited stage (wvi = zvi + hvi)
658    //    wvi- be the h-limited state (wvi- = zvi + hvi-)
659    //
660    //
661    //  where i=0,1,2 denotes the vertex ids
662    //
663    //  Weighted balance between w-limited and h-limited stage is
664    //
665    //    wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi)
666    //
667    //  It follows that the updated wvi is
668    //    wvi := zvi + (1-alpha)*hvi- + alpha*hvi
669    //
670    //   Momentum is balanced between constant and limited
671
672
673    if (alpha < 1) {     
674      for (i=0; i<3; i++) {
675     
676        // FIXME (Ole): Simplify when (if) hvbar gets retired       
677        if (beta_h > epsilon) {   
678          wv[k3+i] = zv[k3+i] + (1-alpha)*hvbar[k3+i] + alpha*hv[i];
679        } else {
680          wv[k3+i] = zv[k3+i] + (1-alpha)*hc_k + alpha*hv[i];   
681        }
682
683        // Update momentum at vertices
684        if (use_centroid_velocities == 1) {
685          // This is a simple, efficient and robust option
686          // It uses first order approximation of velocities, but retains
687          // the order used by stage.
688       
689          // Speeds at centroids
690          if (hc_k > epsilon) {
691            uc = xmomc[k]/hc_k;
692            vc = ymomc[k]/hc_k;
693          } else {
694            uc = 0.0;
695            vc = 0.0;
696          }
697         
698          // Vertex momenta guaranteed to be consistent with depth guaranteeing
699          // controlled speed
700          hv[i] = wv[k3+i] - zv[k3+i]; // Recompute (balanced) vertex depth
701          xmomv[k3+i] = uc*hv[i];
702          ymomv[k3+i] = vc*hv[i];
703         
704        } else {
705          // Update momentum as a linear combination of
706          // xmomc and ymomc (shallow) and momentum
707          // from extrapolator xmomv and ymomv (deep).
708          // This assumes that values from xmomv and ymomv have
709          // been established e.g. by the gradient limiter.
710
711          // FIXME (Ole): I think this should be used with vertex momenta
712          // computed above using centroid_velocities instead of xmomc
713          // and ymomc as they'll be more representative first order
714          // values.
715         
716          xmomv[k3+i] = (1-alpha)*xmomc[k] + alpha*xmomv[k3+i];
717          ymomv[k3+i] = (1-alpha)*ymomc[k] + alpha*ymomv[k3+i];
718       
719        }
720      }
721    }
722  }
723  return 0;
724}
725 
726
727
728
729int _protect(int N,
730             double minimum_allowed_height,
731             double maximum_allowed_speed,
732             double epsilon,
733             double* wc,
734             double* zc,
735             double* xmomc,
736             double* ymomc) {
737
738  int k;
739  double hc;
740  double u, v, reduced_speed;
741
742
743  // Protect against initesimal and negative heights 
744  if (maximum_allowed_speed < epsilon) {
745    for (k=0; k<N; k++) {
746      hc = wc[k] - zc[k];
747
748      if (hc < minimum_allowed_height) {
749       
750        // Set momentum to zero and ensure h is non negative
751        xmomc[k] = 0.0;
752        ymomc[k] = 0.0;
753        if (hc <= 0.0) wc[k] = zc[k];
754      }
755    }
756  } else {
757   
758    // Protect against initesimal and negative heights
759    for (k=0; k<N; k++) {
760      hc = wc[k] - zc[k];
761
762      if (hc < minimum_allowed_height) {
763
764        //New code: Adjust momentum to guarantee speeds are physical
765        //          ensure h is non negative
766             
767        if (hc <= 0.0) {
768                wc[k] = zc[k];
769        xmomc[k] = 0.0;
770        ymomc[k] = 0.0;
771        } else {
772          //Reduce excessive speeds derived from division by small hc
773        //FIXME (Ole): This may be unnecessary with new slope limiters
774        //in effect.
775         
776          u = xmomc[k]/hc;
777          if (fabs(u) > maximum_allowed_speed) {
778            reduced_speed = maximum_allowed_speed * u/fabs(u);
779            //printf("Speed (u) has been reduced from %.3f to %.3f\n",
780            //   u, reduced_speed);
781            xmomc[k] = reduced_speed * hc;
782          }
783
784          v = ymomc[k]/hc;
785          if (fabs(v) > maximum_allowed_speed) {
786            reduced_speed = maximum_allowed_speed * v/fabs(v);
787            //printf("Speed (v) has been reduced from %.3f to %.3f\n",
788            //   v, reduced_speed);
789            ymomc[k] = reduced_speed * hc;
790          }
791        }
792      }
793    }
794  }
795  return 0;
796}
797
798
799
800int _assign_wind_field_values(int N,
801                              double* xmom_update,
802                              double* ymom_update,
803                              double* s_vec,
804                              double* phi_vec,
805                              double cw) {
806
807  // Assign windfield values to momentum updates
808
809  int k;
810  double S, s, phi, u, v;
811
812  for (k=0; k<N; k++) {
813
814    s = s_vec[k];
815    phi = phi_vec[k];
816
817    //Convert to radians
818    phi = phi*pi/180;
819
820    //Compute velocity vector (u, v)
821    u = s*cos(phi);
822    v = s*sin(phi);
823
824    //Compute wind stress
825    S = cw * sqrt(u*u + v*v);
826    xmom_update[k] += S*u;
827    ymom_update[k] += S*v;
828  }
829  return 0;
830}
831
832
833
834///////////////////////////////////////////////////////////////////
835// Gateways to Python
836
837
838
839PyObject *flux_function_central(PyObject *self, PyObject *args) {
840  //
841  // Gateway to innermost flux function.
842  // This is only used by the unit tests as the c implementation is
843  // normally called by compute_fluxes in this module.
844
845
846  PyArrayObject *normal, *ql, *qr,  *edgeflux;
847  double g, epsilon, max_speed, H0, zl, zr;
848
849  if (!PyArg_ParseTuple(args, "OOOddOddd",
850                        &normal, &ql, &qr, &zl, &zr, &edgeflux,
851                        &epsilon, &g, &H0)) {
852    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: flux_function_central could not parse input arguments");
853    return NULL;
854  }
855
856 
857  _flux_function_central((double*) ql -> data, 
858                         (double*) qr -> data, 
859                         zl, 
860                         zr,                                             
861                         ((double*) normal -> data)[0],                                                                         
862                         ((double*) normal -> data)[1],                 
863                         epsilon, H0, g,
864                         (double*) edgeflux -> data, 
865                         &max_speed);
866 
867  return Py_BuildValue("d", max_speed); 
868}
869
870
871
872PyObject *gravity(PyObject *self, PyObject *args) {
873  //
874  //  gravity(g, h, v, x, xmom, ymom)
875  //
876
877
878  PyArrayObject *h, *v, *x, *xmom, *ymom;
879  int k, N, k3, k6;
880  double g, avg_h, zx, zy;
881  double x0, y0, x1, y1, x2, y2, z0, z1, z2;
882  //double epsilon;
883
884  if (!PyArg_ParseTuple(args, "dOOOOO",
885                        &g, &h, &v, &x,
886                        &xmom, &ymom)) {
887    //&epsilon)) {
888    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: gravity could not parse input arguments");
889    return NULL;
890  }
891
892  N = h -> dimensions[0];
893  for (k=0; k<N; k++) {
894    k3 = 3*k;  // base index
895
896    // Get bathymetry
897    z0 = ((double*) v -> data)[k3 + 0];
898    z1 = ((double*) v -> data)[k3 + 1];
899    z2 = ((double*) v -> data)[k3 + 2];
900
901    // Optimise for flat bed
902    // Note (Ole): This didn't produce measurable speed up.
903    // Revisit later
904    //if (fabs(z0-z1)<epsilon && fabs(z1-z2)<epsilon) {
905    //  continue;
906    //}
907
908    // Get average depth from centroid values
909    avg_h = ((double *) h -> data)[k];
910
911    // Compute bed slope
912    k6 = 6*k;  // base index
913 
914    x0 = ((double*) x -> data)[k6 + 0];
915    y0 = ((double*) x -> data)[k6 + 1];
916    x1 = ((double*) x -> data)[k6 + 2];
917    y1 = ((double*) x -> data)[k6 + 3];
918    x2 = ((double*) x -> data)[k6 + 4];
919    y2 = ((double*) x -> data)[k6 + 5];
920
921
922    _gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2, &zx, &zy);
923
924    // Update momentum
925    ((double*) xmom -> data)[k] += -g*zx*avg_h;
926    ((double*) ymom -> data)[k] += -g*zy*avg_h;
927  }
928
929  return Py_BuildValue("");
930}
931
932
933PyObject *manning_friction(PyObject *self, PyObject *args) {
934  //
935  // manning_friction(g, eps, h, uh, vh, eta, xmom_update, ymom_update)
936  //
937
938
939  PyArrayObject *w, *z, *uh, *vh, *eta, *xmom, *ymom;
940  int N;
941  double g, eps;
942
943  if (!PyArg_ParseTuple(args, "ddOOOOOOO",
944                        &g, &eps, &w, &z, &uh, &vh, &eta,
945                        &xmom, &ymom)) {
946    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: manning_friction could not parse input arguments");
947    return NULL;
948  }
949
950
951  N = w -> dimensions[0];
952  _manning_friction(g, eps, N,
953                    (double*) w -> data,
954                    (double*) z -> data,
955                    (double*) uh -> data,
956                    (double*) vh -> data,
957                    (double*) eta -> data,
958                    (double*) xmom -> data,
959                    (double*) ymom -> data);
960
961  return Py_BuildValue("");
962}
963
964
965/*
966PyObject *manning_friction_explicit(PyObject *self, PyObject *args) {
967  //
968  // manning_friction_explicit(g, eps, h, uh, vh, eta, xmom_update, ymom_update)
969  //
970
971
972  PyArrayObject *w, *z, *uh, *vh, *eta, *xmom, *ymom;
973  int N;
974  double g, eps;
975
976  if (!PyArg_ParseTuple(args, "ddOOOOOOO",
977                        &g, &eps, &w, &z, &uh, &vh, &eta,
978                        &xmom, &ymom))
979    return NULL;
980
981  N = w -> dimensions[0];
982  _manning_friction_explicit(g, eps, N,
983                    (double*) w -> data,
984                    (double*) z -> data,
985                    (double*) uh -> data,
986                    (double*) vh -> data,
987                    (double*) eta -> data,
988                    (double*) xmom -> data,
989                    (double*) ymom -> data);
990
991  return Py_BuildValue("");
992}
993*/
994
995
996
997// Computational routine
998int _extrapolate_second_order_sw(int number_of_elements,
999                                  double epsilon,
1000                                  double minimum_allowed_height,
1001                                  double beta_w,
1002                                  double beta_w_dry,
1003                                  double beta_uh,
1004                                  double beta_uh_dry,
1005                                  double beta_vh,
1006                                  double beta_vh_dry,
1007                                  long* surrogate_neighbours,
1008                                  long* number_of_boundaries,
1009                                  double* centroid_coordinates,
1010                                  double* stage_centroid_values,
1011                                  double* xmom_centroid_values,
1012                                  double* ymom_centroid_values,
1013                                  double* elevation_centroid_values,
1014                                  double* vertex_coordinates,
1015                                  double* stage_vertex_values,
1016                                  double* xmom_vertex_values,
1017                                  double* ymom_vertex_values,
1018                                  double* elevation_vertex_values,
1019                                  int optimise_dry_cells) {
1020                                 
1021                                 
1022
1023  // Local variables
1024  double a, b; // Gradient vector used to calculate vertex values from centroids
1025  int k,k0,k1,k2,k3,k6,coord_index,i;
1026  double x,y,x0,y0,x1,y1,x2,y2,xv0,yv0,xv1,yv1,xv2,yv2; // Vertices of the auxiliary triangle
1027  double dx1,dx2,dy1,dy2,dxv0,dxv1,dxv2,dyv0,dyv1,dyv2,dq0,dq1,dq2,area2;
1028  double dqv[3], qmin, qmax, hmin, hmax;
1029  double hc, h0, h1, h2, beta_tmp, hfactor;
1030
1031       
1032  for (k=0; k<number_of_elements; k++) {
1033    k3=k*3;
1034    k6=k*6;
1035
1036   
1037    if (number_of_boundaries[k]==3){
1038      // No neighbours, set gradient on the triangle to zero
1039     
1040      stage_vertex_values[k3]   = stage_centroid_values[k];
1041      stage_vertex_values[k3+1] = stage_centroid_values[k];
1042      stage_vertex_values[k3+2] = stage_centroid_values[k];
1043      xmom_vertex_values[k3]    = xmom_centroid_values[k];
1044      xmom_vertex_values[k3+1]  = xmom_centroid_values[k];
1045      xmom_vertex_values[k3+2]  = xmom_centroid_values[k];
1046      ymom_vertex_values[k3]    = ymom_centroid_values[k];
1047      ymom_vertex_values[k3+1]  = ymom_centroid_values[k];
1048      ymom_vertex_values[k3+2]  = ymom_centroid_values[k];
1049     
1050      continue;
1051    }
1052    else {
1053      // Triangle k has one or more neighbours.
1054      // Get centroid and vertex coordinates of the triangle
1055     
1056      // Get the vertex coordinates
1057      xv0 = vertex_coordinates[k6];   yv0=vertex_coordinates[k6+1];
1058      xv1 = vertex_coordinates[k6+2]; yv1=vertex_coordinates[k6+3];
1059      xv2 = vertex_coordinates[k6+4]; yv2=vertex_coordinates[k6+5];
1060     
1061      // Get the centroid coordinates
1062      coord_index=2*k;
1063      x=centroid_coordinates[coord_index];
1064      y=centroid_coordinates[coord_index+1];
1065     
1066      // Store x- and y- differentials for the vertices of
1067      // triangle k relative to the centroid
1068      dxv0=xv0-x; dxv1=xv1-x; dxv2=xv2-x;
1069      dyv0=yv0-y; dyv1=yv1-y; dyv2=yv2-y;
1070    }
1071
1072
1073           
1074    if (number_of_boundaries[k]<=1){
1075   
1076      //==============================================
1077      // Number of boundaries <= 1
1078      //==============================================   
1079   
1080   
1081      // If no boundaries, auxiliary triangle is formed
1082      // from the centroids of the three neighbours
1083      // If one boundary, auxiliary triangle is formed
1084      // from this centroid and its two neighbours
1085     
1086      k0=surrogate_neighbours[k3];
1087      k1=surrogate_neighbours[k3+1];
1088      k2=surrogate_neighbours[k3+2];
1089     
1090      // Get the auxiliary triangle's vertex coordinates
1091      // (really the centroids of neighbouring triangles)
1092      coord_index=2*k0;
1093      x0=centroid_coordinates[coord_index];
1094      y0=centroid_coordinates[coord_index+1];
1095     
1096      coord_index=2*k1;
1097      x1=centroid_coordinates[coord_index];
1098      y1=centroid_coordinates[coord_index+1];
1099     
1100      coord_index=2*k2;
1101      x2=centroid_coordinates[coord_index];
1102      y2=centroid_coordinates[coord_index+1];
1103     
1104      // Store x- and y- differentials for the vertices
1105      // of the auxiliary triangle
1106      dx1=x1-x0; dx2=x2-x0;
1107      dy1=y1-y0; dy2=y2-y0;
1108     
1109      // Calculate 2*area of the auxiliary triangle
1110      // The triangle is guaranteed to be counter-clockwise     
1111      area2 = dy2*dx1 - dy1*dx2; 
1112     
1113      // If the mesh is 'weird' near the boundary,
1114      // the triangle might be flat or clockwise:
1115      if (area2<=0) {
1116        PyErr_SetString(PyExc_RuntimeError, 
1117                        "shallow_water_ext.c: negative triangle area encountered");
1118        return -1;
1119      } 
1120     
1121      // Calculate heights of neighbouring cells
1122      hc = stage_centroid_values[k]  - elevation_centroid_values[k];
1123      h0 = stage_centroid_values[k0] - elevation_centroid_values[k0];
1124      h1 = stage_centroid_values[k1] - elevation_centroid_values[k1];
1125      h2 = stage_centroid_values[k2] - elevation_centroid_values[k2];
1126      hmin = min(min(h0,min(h1,h2)),hc);
1127      //hfactor = hc/(hc + 1.0);
1128
1129      hfactor = 0.0;
1130      if (hmin > 0.001 ) {
1131          hfactor = (hmin-0.001)/(hmin+0.004);
1132      }
1133     
1134      if (optimise_dry_cells) {     
1135        // Check if linear reconstruction is necessary for triangle k
1136        // This check will exclude dry cells.
1137
1138        hmax = max(h0,max(h1,h2));     
1139        if (hmax < epsilon) {
1140          continue;
1141        }
1142      }
1143
1144           
1145      //-----------------------------------
1146      // stage
1147      //-----------------------------------     
1148     
1149      // Calculate the difference between vertex 0 of the auxiliary
1150      // triangle and the centroid of triangle k
1151      dq0=stage_centroid_values[k0]-stage_centroid_values[k];
1152     
1153      // Calculate differentials between the vertices
1154      // of the auxiliary triangle (centroids of neighbouring triangles)
1155      dq1=stage_centroid_values[k1]-stage_centroid_values[k0];
1156      dq2=stage_centroid_values[k2]-stage_centroid_values[k0];
1157     
1158      // Calculate the gradient of stage on the auxiliary triangle
1159      a = dy2*dq1 - dy1*dq2;
1160      a /= area2;
1161      b = dx1*dq2 - dx2*dq1;
1162      b /= area2;
1163     
1164      // Calculate provisional jumps in stage from the centroid
1165      // of triangle k to its vertices, to be limited
1166      dqv[0]=a*dxv0+b*dyv0;
1167      dqv[1]=a*dxv1+b*dyv1;
1168      dqv[2]=a*dxv2+b*dyv2;
1169     
1170      // Now we want to find min and max of the centroid and the
1171      // vertices of the auxiliary triangle and compute jumps
1172      // from the centroid to the min and max
1173      find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
1174     
1175      // Playing with dry wet interface
1176      //hmin = qmin;
1177      //beta_tmp = beta_w_dry;
1178      //if (hmin>minimum_allowed_height)
1179      beta_tmp = beta_w_dry + (beta_w - beta_w_dry) * hfactor;
1180       
1181      //printf("min_alled_height = %f\n",minimum_allowed_height);
1182      //printf("hmin = %f\n",hmin);
1183      //printf("beta_w = %f\n",beta_w);
1184      //printf("beta_tmp = %f\n",beta_tmp);
1185      // Limit the gradient
1186      limit_gradient(dqv,qmin,qmax,beta_tmp); 
1187     
1188      for (i=0;i<3;i++)
1189        stage_vertex_values[k3+i]=stage_centroid_values[k]+dqv[i];
1190     
1191     
1192      //-----------------------------------
1193      // xmomentum
1194      //-----------------------------------           
1195
1196      // Calculate the difference between vertex 0 of the auxiliary
1197      // triangle and the centroid of triangle k     
1198      dq0=xmom_centroid_values[k0]-xmom_centroid_values[k];
1199     
1200      // Calculate differentials between the vertices
1201      // of the auxiliary triangle
1202      dq1=xmom_centroid_values[k1]-xmom_centroid_values[k0];
1203      dq2=xmom_centroid_values[k2]-xmom_centroid_values[k0];
1204     
1205      // Calculate the gradient of xmom on the auxiliary triangle
1206      a = dy2*dq1 - dy1*dq2;
1207      a /= area2;
1208      b = dx1*dq2 - dx2*dq1;
1209      b /= area2;
1210     
1211      // Calculate provisional jumps in stage from the centroid
1212      // of triangle k to its vertices, to be limited     
1213      dqv[0]=a*dxv0+b*dyv0;
1214      dqv[1]=a*dxv1+b*dyv1;
1215      dqv[2]=a*dxv2+b*dyv2;
1216     
1217      // Now we want to find min and max of the centroid and the
1218      // vertices of the auxiliary triangle and compute jumps
1219      // from the centroid to the min and max
1220      find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
1221      //beta_tmp = beta_uh;
1222      //if (hmin<minimum_allowed_height)
1223      //beta_tmp = beta_uh_dry;
1224      beta_tmp = beta_uh_dry + (beta_uh - beta_uh_dry) * hfactor;
1225
1226      // Limit the gradient
1227      limit_gradient(dqv,qmin,qmax,beta_tmp);
1228
1229      for (i=0;i<3;i++)
1230        xmom_vertex_values[k3+i]=xmom_centroid_values[k]+dqv[i];
1231     
1232     
1233      //-----------------------------------
1234      // ymomentum
1235      //-----------------------------------                 
1236
1237      // Calculate the difference between vertex 0 of the auxiliary
1238      // triangle and the centroid of triangle k
1239      dq0=ymom_centroid_values[k0]-ymom_centroid_values[k];
1240     
1241      // Calculate differentials between the vertices
1242      // of the auxiliary triangle
1243      dq1=ymom_centroid_values[k1]-ymom_centroid_values[k0];
1244      dq2=ymom_centroid_values[k2]-ymom_centroid_values[k0];
1245     
1246      // Calculate the gradient of xmom on the auxiliary triangle
1247      a = dy2*dq1 - dy1*dq2;
1248      a /= area2;
1249      b = dx1*dq2 - dx2*dq1;
1250      b /= area2;
1251     
1252      // Calculate provisional jumps in stage from the centroid
1253      // of triangle k to its vertices, to be limited
1254      dqv[0]=a*dxv0+b*dyv0;
1255      dqv[1]=a*dxv1+b*dyv1;
1256      dqv[2]=a*dxv2+b*dyv2;
1257     
1258      // Now we want to find min and max of the centroid and the
1259      // vertices of the auxiliary triangle and compute jumps
1260      // from the centroid to the min and max
1261      find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
1262     
1263      //beta_tmp = beta_vh;
1264      //
1265      //if (hmin<minimum_allowed_height)
1266      //beta_tmp = beta_vh_dry;
1267      beta_tmp = beta_vh_dry + (beta_vh - beta_vh_dry) * hfactor;       
1268
1269      // Limit the gradient
1270      limit_gradient(dqv,qmin,qmax,beta_tmp);
1271     
1272      for (i=0;i<3;i++)
1273        ymom_vertex_values[k3+i]=ymom_centroid_values[k]+dqv[i];
1274       
1275    } // End number_of_boundaries <=1
1276    else{
1277
1278      //==============================================
1279      // Number of boundaries == 2
1280      //==============================================       
1281       
1282      // One internal neighbour and gradient is in direction of the neighbour's centroid
1283     
1284      // Find the only internal neighbour (k1?)
1285      for (k2=k3;k2<k3+3;k2++){
1286        // Find internal neighbour of triangle k     
1287        // k2 indexes the edges of triangle k   
1288     
1289        if (surrogate_neighbours[k2]!=k)
1290          break;
1291      }
1292     
1293      if ((k2==k3+3)) {
1294        // If we didn't find an internal neighbour
1295        PyErr_SetString(PyExc_RuntimeError, 
1296                        "shallow_water_ext.c: Internal neighbour not found");     
1297        return -1;
1298      }
1299     
1300      k1=surrogate_neighbours[k2];
1301     
1302      // The coordinates of the triangle are already (x,y).
1303      // Get centroid of the neighbour (x1,y1)
1304      coord_index=2*k1;
1305      x1=centroid_coordinates[coord_index];
1306      y1=centroid_coordinates[coord_index+1];
1307     
1308      // Compute x- and y- distances between the centroid of
1309      // triangle k and that of its neighbour
1310      dx1=x1-x; dy1=y1-y;
1311     
1312      // Set area2 as the square of the distance
1313      area2=dx1*dx1+dy1*dy1;
1314     
1315      // Set dx2=(x1-x0)/((x1-x0)^2+(y1-y0)^2)
1316      // and dy2=(y1-y0)/((x1-x0)^2+(y1-y0)^2) which
1317      // respectively correspond to the x- and y- gradients
1318      // of the conserved quantities
1319      dx2=1.0/area2;
1320      dy2=dx2*dy1;
1321      dx2*=dx1;
1322     
1323     
1324      //-----------------------------------
1325      // stage
1326      //-----------------------------------           
1327
1328      // Compute differentials
1329      dq1=stage_centroid_values[k1]-stage_centroid_values[k];
1330     
1331      // Calculate the gradient between the centroid of triangle k
1332      // and that of its neighbour
1333      a=dq1*dx2;
1334      b=dq1*dy2;
1335     
1336      // Calculate provisional vertex jumps, to be limited
1337      dqv[0]=a*dxv0+b*dyv0;
1338      dqv[1]=a*dxv1+b*dyv1;
1339      dqv[2]=a*dxv2+b*dyv2;
1340     
1341      // Now limit the jumps
1342      if (dq1>=0.0){
1343        qmin=0.0;
1344        qmax=dq1;
1345      }
1346      else{
1347        qmin=dq1;
1348        qmax=0.0;
1349      }
1350     
1351      // Limit the gradient
1352      limit_gradient(dqv,qmin,qmax,beta_w);
1353     
1354      for (i=0;i<3;i++)
1355        stage_vertex_values[k3+i]=stage_centroid_values[k]+dqv[i];
1356     
1357     
1358      //-----------------------------------
1359      // xmomentum
1360      //-----------------------------------                       
1361     
1362      // Compute differentials
1363      dq1=xmom_centroid_values[k1]-xmom_centroid_values[k];
1364     
1365      // Calculate the gradient between the centroid of triangle k
1366      // and that of its neighbour
1367      a=dq1*dx2;
1368      b=dq1*dy2;
1369     
1370      // Calculate provisional vertex jumps, to be limited
1371      dqv[0]=a*dxv0+b*dyv0;
1372      dqv[1]=a*dxv1+b*dyv1;
1373      dqv[2]=a*dxv2+b*dyv2;
1374     
1375      // Now limit the jumps
1376      if (dq1>=0.0){
1377        qmin=0.0;
1378        qmax=dq1;
1379      }
1380      else{
1381        qmin=dq1;
1382        qmax=0.0;
1383      }
1384     
1385      // Limit the gradient     
1386      limit_gradient(dqv,qmin,qmax,beta_w);
1387     
1388      for (i=0;i<3;i++)
1389        xmom_vertex_values[k3+i]=xmom_centroid_values[k]+dqv[i];
1390     
1391     
1392      //-----------------------------------
1393      // ymomentum
1394      //-----------------------------------                       
1395
1396      // Compute differentials
1397      dq1=ymom_centroid_values[k1]-ymom_centroid_values[k];
1398     
1399      // Calculate the gradient between the centroid of triangle k
1400      // and that of its neighbour
1401      a=dq1*dx2;
1402      b=dq1*dy2;
1403     
1404      // Calculate provisional vertex jumps, to be limited
1405      dqv[0]=a*dxv0+b*dyv0;
1406      dqv[1]=a*dxv1+b*dyv1;
1407      dqv[2]=a*dxv2+b*dyv2;
1408     
1409      // Now limit the jumps
1410      if (dq1>=0.0){
1411        qmin=0.0;
1412        qmax=dq1;
1413      }
1414      else{
1415        qmin=dq1;
1416        qmax=0.0;
1417      }
1418     
1419      // Limit the gradient
1420      limit_gradient(dqv,qmin,qmax,beta_w);
1421     
1422      for (i=0;i<3;i++)
1423        ymom_vertex_values[k3+i]=ymom_centroid_values[k]+dqv[i];
1424       
1425    } // else [number_of_boundaries==2]
1426  } // for k=0 to number_of_elements-1
1427 
1428  return 0;
1429}                       
1430
1431
1432PyObject *extrapolate_second_order_sw(PyObject *self, PyObject *args) {
1433  /*Compute the vertex values based on a linear reconstruction
1434    on each triangle
1435   
1436    These values are calculated as follows:
1437    1) For each triangle not adjacent to a boundary, we consider the
1438       auxiliary triangle formed by the centroids of its three
1439       neighbours.
1440    2) For each conserved quantity, we integrate around the auxiliary
1441       triangle's boundary the product of the quantity and the outward
1442       normal vector. Dividing by the triangle area gives (a,b), the
1443       average of the vector (q_x,q_y) on the auxiliary triangle.
1444       We suppose that the linear reconstruction on the original
1445       triangle has gradient (a,b).
1446    3) Provisional vertex jumps dqv[0,1,2] are computed and these are
1447       then limited by calling the functions find_qmin_and_qmax and
1448       limit_gradient
1449
1450    Python call:
1451    extrapolate_second_order_sw(domain.surrogate_neighbours,
1452                                domain.number_of_boundaries
1453                                domain.centroid_coordinates,
1454                                Stage.centroid_values
1455                                Xmom.centroid_values
1456                                Ymom.centroid_values
1457                                domain.vertex_coordinates,
1458                                Stage.vertex_values,
1459                                Xmom.vertex_values,
1460                                Ymom.vertex_values)
1461
1462    Post conditions:
1463            The vertices of each triangle have values from a
1464            limited linear reconstruction
1465            based on centroid values
1466
1467  */
1468  PyArrayObject *surrogate_neighbours,
1469    *number_of_boundaries,
1470    *centroid_coordinates,
1471    *stage_centroid_values,
1472    *xmom_centroid_values,
1473    *ymom_centroid_values,
1474    *elevation_centroid_values,
1475    *vertex_coordinates,
1476    *stage_vertex_values,
1477    *xmom_vertex_values,
1478    *ymom_vertex_values,
1479    *elevation_vertex_values;
1480 
1481  PyObject *domain;
1482
1483 
1484  double beta_w, beta_w_dry, beta_uh, beta_uh_dry, beta_vh, beta_vh_dry;   
1485  double minimum_allowed_height, epsilon;
1486  int optimise_dry_cells, number_of_elements, e;
1487 
1488  // Provisional jumps from centroids to v'tices and safety factor re limiting
1489  // by which these jumps are limited
1490  // Convert Python arguments to C
1491  if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOi",
1492                        &domain,
1493                        &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                        &optimise_dry_cells)) {                 
1506                       
1507    PyErr_SetString(PyExc_RuntimeError, 
1508                    "Input arguments to extrapolate_second_order_sw failed");
1509    return NULL;
1510  }
1511
1512  // Get the safety factor beta_w, set in the config.py file.
1513  // This is used in the limiting process
1514 
1515
1516  beta_w                 = get_python_double(domain,"beta_w");
1517  beta_w_dry             = get_python_double(domain,"beta_w_dry");
1518  beta_uh                = get_python_double(domain,"beta_uh");
1519  beta_uh_dry            = get_python_double(domain,"beta_uh_dry");
1520  beta_vh                = get_python_double(domain,"beta_vh");
1521  beta_vh_dry            = get_python_double(domain,"beta_vh_dry"); 
1522
1523  minimum_allowed_height = get_python_double(domain,"minimum_allowed_height");
1524  epsilon                = get_python_double(domain,"epsilon");
1525
1526  number_of_elements = stage_centroid_values -> dimensions[0]; 
1527
1528  // Call underlying computational routine
1529  e = _extrapolate_second_order_sw(number_of_elements,
1530                                   epsilon,
1531                                   minimum_allowed_height,
1532                                   beta_w,
1533                                   beta_w_dry,
1534                                   beta_uh,
1535                                   beta_uh_dry,
1536                                   beta_vh,
1537                                   beta_vh_dry,
1538                                   (long*) surrogate_neighbours -> data,
1539                                   (long*) number_of_boundaries -> data,
1540                                   (double*) centroid_coordinates -> data,
1541                                   (double*) stage_centroid_values -> data,
1542                                   (double*) xmom_centroid_values -> data,
1543                                   (double*) ymom_centroid_values -> data,
1544                                   (double*) elevation_centroid_values -> data,
1545                                   (double*) vertex_coordinates -> data,
1546                                   (double*) stage_vertex_values -> data,
1547                                   (double*) xmom_vertex_values -> data,
1548                                   (double*) ymom_vertex_values -> data,
1549                                   (double*) elevation_vertex_values -> data,
1550                                   optimise_dry_cells);
1551  if (e == -1) {
1552    // Use error string set inside computational routine
1553    return NULL;
1554  }                               
1555 
1556 
1557  return Py_BuildValue("");
1558 
1559}// extrapolate_second-order_sw
1560
1561
1562
1563
1564PyObject *rotate(PyObject *self, PyObject *args, PyObject *kwargs) {
1565  //
1566  // r = rotate(q, normal, direction=1)
1567  //
1568  // Where q is assumed to be a Float numeric array of length 3 and
1569  // normal a Float numeric array of length 2.
1570
1571  // FIXME(Ole): I don't think this is used anymore
1572
1573  PyObject *Q, *Normal;
1574  PyArrayObject *q, *r, *normal;
1575
1576  static char *argnames[] = {"q", "normal", "direction", NULL};
1577  int dimensions[1], i, direction=1;
1578  double n1, n2;
1579
1580  // Convert Python arguments to C
1581  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OO|i", argnames,
1582                                   &Q, &Normal, &direction)) {
1583    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: rotate could not parse input arguments");
1584    return NULL;
1585  } 
1586
1587  // Input checks (convert sequences into numeric arrays)
1588  q = (PyArrayObject *)
1589    PyArray_ContiguousFromObject(Q, PyArray_DOUBLE, 0, 0);
1590  normal = (PyArrayObject *)
1591    PyArray_ContiguousFromObject(Normal, PyArray_DOUBLE, 0, 0);
1592
1593
1594  if (normal -> dimensions[0] != 2) {
1595    PyErr_SetString(PyExc_RuntimeError, "Normal vector must have 2 components");
1596    return NULL;
1597  }
1598
1599  // Allocate space for return vector r (don't DECREF)
1600  dimensions[0] = 3;
1601  r = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
1602
1603  // Copy
1604  for (i=0; i<3; i++) {
1605    ((double *) (r -> data))[i] = ((double *) (q -> data))[i];
1606  }
1607
1608  // Get normal and direction
1609  n1 = ((double *) normal -> data)[0];
1610  n2 = ((double *) normal -> data)[1];
1611  if (direction == -1) n2 = -n2;
1612
1613  // Rotate
1614  _rotate((double *) r -> data, n1, n2);
1615
1616  // Release numeric arrays
1617  Py_DECREF(q);
1618  Py_DECREF(normal);
1619
1620  // Return result using PyArray to avoid memory leak
1621  return PyArray_Return(r);
1622}
1623
1624
1625// Computational function for flux computation
1626double _compute_fluxes_central(int number_of_elements,
1627                               double timestep,
1628                               double epsilon,
1629                               double H0,
1630                               double g,
1631                               long* neighbours,
1632                               long* neighbour_edges,
1633                               double* normals,
1634                               double* edgelengths, 
1635                               double* radii, 
1636                               double* areas,
1637                               long* tri_full_flag,
1638                               double* stage_edge_values,
1639                               double* xmom_edge_values,
1640                               double* ymom_edge_values,
1641                               double* bed_edge_values,
1642                               double* stage_boundary_values,
1643                               double* xmom_boundary_values,
1644                               double* ymom_boundary_values,
1645                               double* stage_explicit_update,
1646                               double* xmom_explicit_update,
1647                               double* ymom_explicit_update,
1648                               long* already_computed_flux,
1649                               double* max_speed_array,
1650                               int optimise_dry_cells) {
1651                               
1652  // Local variables
1653  double max_speed, length, area, zl, zr;
1654  int k, i, m, n;
1655  int ki, nm=0, ki2; // Index shorthands
1656 
1657  // Workspace (making them static actually made function slightly slower (Ole)) 
1658  double ql[3], qr[3], edgeflux[3]; // Work array for summing up fluxes
1659
1660  static long call=1; // Static local variable flagging already computed flux
1661                               
1662       
1663  // Start computation
1664  call++; // Flag 'id' of flux calculation for this timestep
1665 
1666  // Set explicit_update to zero for all conserved_quantities.
1667  // This assumes compute_fluxes called before forcing terms
1668  for (k=0; k<number_of_elements; k++) {
1669    stage_explicit_update[k]=0.0;
1670    xmom_explicit_update[k]=0.0;
1671    ymom_explicit_update[k]=0.0; 
1672  }
1673
1674  // For all triangles
1675  for (k=0; k<number_of_elements; k++) {
1676   
1677    // Loop through neighbours and compute edge flux for each 
1678    for (i=0; i<3; i++) {
1679      ki = k*3+i; // Linear index (triangle k, edge i)
1680     
1681      if (already_computed_flux[ki] == call)
1682        // We've already computed the flux across this edge
1683        continue;
1684       
1685       
1686      ql[0] = stage_edge_values[ki];
1687      ql[1] = xmom_edge_values[ki];
1688      ql[2] = ymom_edge_values[ki];
1689      zl = bed_edge_values[ki];
1690
1691      // Quantities at neighbour on nearest face
1692      n = neighbours[ki];
1693      if (n < 0) {
1694        // Neighbour is a boundary condition
1695        m = -n-1; // Convert negative flag to boundary index
1696       
1697        qr[0] = stage_boundary_values[m];
1698        qr[1] = xmom_boundary_values[m];
1699        qr[2] = ymom_boundary_values[m];
1700        zr = zl; // Extend bed elevation to boundary
1701      } else {
1702        // Neighbour is a real element
1703        m = neighbour_edges[ki];
1704        nm = n*3+m; // Linear index (triangle n, edge m)
1705       
1706        qr[0] = stage_edge_values[nm];
1707        qr[1] = xmom_edge_values[nm];
1708        qr[2] = ymom_edge_values[nm];
1709        zr = bed_edge_values[nm];
1710      }
1711     
1712     
1713      if (optimise_dry_cells) {     
1714        // Check if flux calculation is necessary across this edge
1715        // This check will exclude dry cells.
1716        // This will also optimise cases where zl != zr as
1717        // long as both are dry
1718
1719        if ( fabs(ql[0] - zl) < epsilon && 
1720             fabs(qr[0] - zr) < epsilon ) {
1721          // Cell boundary is dry
1722         
1723          already_computed_flux[ki] = call; // #k Done 
1724          if (n>=0)
1725            already_computed_flux[nm] = call; // #n Done
1726       
1727          max_speed = 0.0;
1728          continue;
1729        }
1730      }
1731   
1732     
1733      // Outward pointing normal vector (domain.normals[k, 2*i:2*i+2])
1734      ki2 = 2*ki; //k*6 + i*2
1735
1736      // Edge flux computation (triangle k, edge i)
1737      _flux_function_central(ql, qr, zl, zr,
1738                             normals[ki2], normals[ki2+1],
1739                             epsilon, H0, g,
1740                             edgeflux, &max_speed);
1741     
1742     
1743      // Multiply edgeflux by edgelength
1744      length = edgelengths[ki];
1745      edgeflux[0] *= length;           
1746      edgeflux[1] *= length;           
1747      edgeflux[2] *= length;                       
1748     
1749     
1750      // Update triangle k with flux from edge i
1751      stage_explicit_update[k] -= edgeflux[0];
1752      xmom_explicit_update[k] -= edgeflux[1];
1753      ymom_explicit_update[k] -= edgeflux[2];
1754     
1755      already_computed_flux[ki] = call; // #k Done
1756     
1757     
1758      // Update neighbour n with same flux but reversed sign
1759      if (n>=0) {
1760        stage_explicit_update[n] += edgeflux[0];
1761        xmom_explicit_update[n] += edgeflux[1];
1762        ymom_explicit_update[n] += edgeflux[2];
1763       
1764        already_computed_flux[nm] = call; // #n Done
1765      }
1766
1767
1768      // Update timestep based on edge i and possibly neighbour n
1769      if (tri_full_flag[k] == 1) {
1770        if (max_speed > epsilon) {
1771       
1772          // Original CFL calculation
1773          timestep = min(timestep, radii[k]/max_speed);
1774          if (n>=0) 
1775            timestep = min(timestep, radii[n]/max_speed);
1776         
1777          // Ted Rigby's suggested less conservative version
1778          //if (n>=0) {             
1779          //  timestep = min(timestep, (radii[k]+radii[n])/max_speed);
1780          //} else {
1781          //  timestep = min(timestep, radii[k]/max_speed);
1782          // }
1783        }
1784      }
1785     
1786    } // End edge i (and neighbour n)
1787   
1788   
1789    // Normalise triangle k by area and store for when all conserved
1790    // quantities get updated
1791    area = areas[k];
1792    stage_explicit_update[k] /= area;
1793    xmom_explicit_update[k] /= area;
1794    ymom_explicit_update[k] /= area;
1795   
1796   
1797    // Keep track of maximal speeds
1798    max_speed_array[k] = max_speed;
1799   
1800  } // End triangle k
1801       
1802                               
1803                               
1804  return timestep;
1805}
1806
1807//=========================================================================
1808// Python Glue
1809//=========================================================================
1810
1811PyObject *compute_fluxes_ext_central_new(PyObject *self, PyObject *args) {
1812  /*Compute all fluxes and the timestep suitable for all volumes
1813    in domain.
1814
1815    Compute total flux for each conserved quantity using "flux_function_central"
1816
1817    Fluxes across each edge are scaled by edgelengths and summed up
1818    Resulting flux is then scaled by area and stored in
1819    explicit_update for each of the three conserved quantities
1820    stage, xmomentum and ymomentum
1821
1822    The maximal allowable speed computed by the flux_function for each volume
1823    is converted to a timestep that must not be exceeded. The minimum of
1824    those is computed as the next overall timestep.
1825
1826    Python call:
1827    timestep = compute_fluxes(timestep, domain, stage, xmom, ymom, bed)
1828
1829
1830    Post conditions:
1831      domain.explicit_update is reset to computed flux values
1832      returns timestep which is the largest step satisfying all volumes.
1833
1834
1835  */
1836
1837    PyObject
1838        *domain,
1839        *stage, 
1840        *xmom, 
1841        *ymom, 
1842        *bed;
1843
1844    PyArrayObject
1845        *neighbours, 
1846        *neighbour_edges,
1847        *normals, 
1848        *edgelengths, 
1849        *radii, 
1850        *areas,
1851        *tri_full_flag,
1852        *stage_edge_values,
1853        *xmom_edge_values,
1854        *ymom_edge_values,
1855        *bed_edge_values,
1856        *stage_boundary_values,
1857        *xmom_boundary_values,
1858        *ymom_boundary_values,
1859        *stage_explicit_update,
1860        *xmom_explicit_update,
1861        *ymom_explicit_update,
1862        *already_computed_flux, //Tracks whether the flux across an edge has already been computed
1863        *max_speed_array; //Keeps track of max speeds for each triangle
1864
1865   
1866    double timestep, epsilon, H0, g;
1867    int optimise_dry_cells;
1868   
1869    // Convert Python arguments to C
1870    if (!PyArg_ParseTuple(args, "dOOOO", &timestep, &domain, &stage, &xmom, &ymom, &bed )) {
1871        PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
1872        return NULL;
1873    }
1874
1875    epsilon           = get_python_double(domain,"epsilon");
1876    H0                = get_python_double(domain,"H0");
1877    g                 = get_python_double(domain,"g");
1878    optimise_dry_cells = get_python_integer(domain,"optimse_dry_cells");
1879   
1880    neighbours        = get_consecutive_array(domain, "neighbours");
1881    neighbour_edges   = get_consecutive_array(domain, "neighbour_edges"); 
1882    normals           = get_consecutive_array(domain, "normals");
1883    edgelengths       = get_consecutive_array(domain, "edge_lengths");   
1884    radii             = get_consecutive_array(domain, "radii");   
1885    areas             = get_consecutive_array(domain, "areas");   
1886    tri_full_flag     = get_consecutive_array(domain, "normals");
1887    already_computed_flux  = get_consecutive_array(domain, "already_computed_flux");
1888    max_speed_array   = get_consecutive_array(domain, "max_speed");
1889   
1890    stage_edge_values = get_consecutive_array(stage, "edge_values");   
1891    xmom_edge_values  = get_consecutive_array(xmom, "edge_values");   
1892    ymom_edge_values  = get_consecutive_array(ymom, "edge_values"); 
1893    bed_edge_values   = get_consecutive_array(bed, "edge_values");   
1894
1895    stage_boundary_values = get_consecutive_array(stage, "boundary_values");   
1896    xmom_boundary_values  = get_consecutive_array(xmom, "boundary_values");   
1897    ymom_boundary_values  = get_consecutive_array(ymom, "boundary_values"); 
1898
1899    stage_explicit_update = get_consecutive_array(stage, "explicit_update");   
1900    xmom_explicit_update  = get_consecutive_array(xmom, "explicit_update");   
1901    ymom_explicit_update  = get_consecutive_array(ymom, "explicit_update"); 
1902
1903
1904  int number_of_elements = stage_edge_values -> dimensions[0];
1905
1906  // Call underlying flux computation routine and update
1907  // the explicit update arrays
1908  timestep = _compute_fluxes_central(number_of_elements,
1909                                     timestep,
1910                                     epsilon,
1911                                     H0,
1912                                     g,
1913                                     (long*) neighbours -> data,
1914                                     (long*) neighbour_edges -> data,
1915                                     (double*) normals -> data,
1916                                     (double*) edgelengths -> data, 
1917                                     (double*) radii -> data, 
1918                                     (double*) areas -> data,
1919                                     (long*) tri_full_flag -> data,
1920                                     (double*) stage_edge_values -> data,
1921                                     (double*) xmom_edge_values -> data,
1922                                     (double*) ymom_edge_values -> data,
1923                                     (double*) bed_edge_values -> data,
1924                                     (double*) stage_boundary_values -> data,
1925                                     (double*) xmom_boundary_values -> data,
1926                                     (double*) ymom_boundary_values -> data,
1927                                     (double*) stage_explicit_update -> data,
1928                                     (double*) xmom_explicit_update -> data,
1929                                     (double*) ymom_explicit_update -> data,
1930                                     (long*) already_computed_flux -> data,
1931                                     (double*) max_speed_array -> data,
1932                                     optimise_dry_cells);
1933
1934  Py_DECREF(neighbours);
1935  Py_DECREF(neighbour_edges);
1936  Py_DECREF(normals);
1937  Py_DECREF(edgelengths);
1938  Py_DECREF(radii);
1939  Py_DECREF(areas);
1940  Py_DECREF(tri_full_flag);
1941  Py_DECREF(already_computed_flux);
1942  Py_DECREF(max_speed_array);
1943  Py_DECREF(stage_edge_values);
1944  Py_DECREF(xmom_edge_values);
1945  Py_DECREF(ymom_edge_values);
1946  Py_DECREF(bed_edge_values);
1947  Py_DECREF(stage_boundary_values);
1948  Py_DECREF(xmom_boundary_values);
1949  Py_DECREF(ymom_boundary_values);
1950  Py_DECREF(stage_explicit_update);
1951  Py_DECREF(xmom_explicit_update);
1952  Py_DECREF(ymom_explicit_update);
1953
1954 
1955  // Return updated flux timestep
1956  return Py_BuildValue("d", timestep);
1957}
1958
1959
1960
1961
1962
1963
1964PyObject *compute_fluxes_ext_central(PyObject *self, PyObject *args) {
1965  /*Compute all fluxes and the timestep suitable for all volumes
1966    in domain.
1967
1968    Compute total flux for each conserved quantity using "flux_function_central"
1969
1970    Fluxes across each edge are scaled by edgelengths and summed up
1971    Resulting flux is then scaled by area and stored in
1972    explicit_update for each of the three conserved quantities
1973    stage, xmomentum and ymomentum
1974
1975    The maximal allowable speed computed by the flux_function for each volume
1976    is converted to a timestep that must not be exceeded. The minimum of
1977    those is computed as the next overall timestep.
1978
1979    Python call:
1980    domain.timestep = compute_fluxes(timestep,
1981                                     domain.epsilon,
1982                                     domain.H0,
1983                                     domain.g,
1984                                     domain.neighbours,
1985                                     domain.neighbour_edges,
1986                                     domain.normals,
1987                                     domain.edgelengths,
1988                                     domain.radii,
1989                                     domain.areas,
1990                                     tri_full_flag,
1991                                     Stage.edge_values,
1992                                     Xmom.edge_values,
1993                                     Ymom.edge_values,
1994                                     Bed.edge_values,
1995                                     Stage.boundary_values,
1996                                     Xmom.boundary_values,
1997                                     Ymom.boundary_values,
1998                                     Stage.explicit_update,
1999                                     Xmom.explicit_update,
2000                                     Ymom.explicit_update,
2001                                     already_computed_flux,
2002                                     optimise_dry_cells)                                     
2003
2004
2005    Post conditions:
2006      domain.explicit_update is reset to computed flux values
2007      domain.timestep is set to the largest step satisfying all volumes.
2008
2009
2010  */
2011
2012
2013  PyArrayObject *neighbours, *neighbour_edges,
2014    *normals, *edgelengths, *radii, *areas,
2015    *tri_full_flag,
2016    *stage_edge_values,
2017    *xmom_edge_values,
2018    *ymom_edge_values,
2019    *bed_edge_values,
2020    *stage_boundary_values,
2021    *xmom_boundary_values,
2022    *ymom_boundary_values,
2023    *stage_explicit_update,
2024    *xmom_explicit_update,
2025    *ymom_explicit_update,
2026    *already_computed_flux, //Tracks whether the flux across an edge has already been computed
2027    *max_speed_array; //Keeps track of max speeds for each triangle
2028
2029   
2030  double timestep, epsilon, H0, g;
2031  int optimise_dry_cells;
2032   
2033  // Convert Python arguments to C
2034  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOi",
2035                        &timestep,
2036                        &epsilon,
2037                        &H0,
2038                        &g,
2039                        &neighbours,
2040                        &neighbour_edges,
2041                        &normals,
2042                        &edgelengths, &radii, &areas,
2043                        &tri_full_flag,
2044                        &stage_edge_values,
2045                        &xmom_edge_values,
2046                        &ymom_edge_values,
2047                        &bed_edge_values,
2048                        &stage_boundary_values,
2049                        &xmom_boundary_values,
2050                        &ymom_boundary_values,
2051                        &stage_explicit_update,
2052                        &xmom_explicit_update,
2053                        &ymom_explicit_update,
2054                        &already_computed_flux,
2055                        &max_speed_array,
2056                        &optimise_dry_cells)) {
2057    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
2058    return NULL;
2059  }
2060
2061 
2062  int number_of_elements = stage_edge_values -> dimensions[0];
2063
2064  // Call underlying flux computation routine and update
2065  // the explicit update arrays
2066  timestep = _compute_fluxes_central(number_of_elements,
2067                                     timestep,
2068                                     epsilon,
2069                                     H0,
2070                                     g,
2071                                     (long*) neighbours -> data,
2072                                     (long*) neighbour_edges -> data,
2073                                     (double*) normals -> data,
2074                                     (double*) edgelengths -> data, 
2075                                     (double*) radii -> data, 
2076                                     (double*) areas -> data,
2077                                     (long*) tri_full_flag -> data,
2078                                     (double*) stage_edge_values -> data,
2079                                     (double*) xmom_edge_values -> data,
2080                                     (double*) ymom_edge_values -> data,
2081                                     (double*) bed_edge_values -> data,
2082                                     (double*) stage_boundary_values -> data,
2083                                     (double*) xmom_boundary_values -> data,
2084                                     (double*) ymom_boundary_values -> data,
2085                                     (double*) stage_explicit_update -> data,
2086                                     (double*) xmom_explicit_update -> data,
2087                                     (double*) ymom_explicit_update -> data,
2088                                     (long*) already_computed_flux -> data,
2089                                     (double*) max_speed_array -> data,
2090                                     optimise_dry_cells);
2091 
2092  // Return updated flux timestep
2093  return Py_BuildValue("d", timestep);
2094}
2095
2096
2097
2098
2099
2100PyObject *compute_fluxes_ext_kinetic(PyObject *self, PyObject *args) {
2101  /*Compute all fluxes and the timestep suitable for all volumes
2102    in domain.
2103
2104    THIS IS AN EXPERIMENTAL FUNCTION - NORMALLY flux_function_central IS USED.
2105   
2106    Compute total flux for each conserved quantity using "flux_function_kinetic"
2107
2108    Fluxes across each edge are scaled by edgelengths and summed up
2109    Resulting flux is then scaled by area and stored in
2110    explicit_update for each of the three conserved quantities
2111    stage, xmomentum and ymomentum
2112
2113    The maximal allowable speed computed by the flux_function for each volume
2114    is converted to a timestep that must not be exceeded. The minimum of
2115    those is computed as the next overall timestep.
2116
2117    Python call:
2118    domain.timestep = compute_fluxes(timestep,
2119                                     domain.epsilon,
2120                                     domain.H0,
2121                                     domain.g,
2122                                     domain.neighbours,
2123                                     domain.neighbour_edges,
2124                                     domain.normals,
2125                                     domain.edgelengths,
2126                                     domain.radii,
2127                                     domain.areas,
2128                                     Stage.edge_values,
2129                                     Xmom.edge_values,
2130                                     Ymom.edge_values,
2131                                     Bed.edge_values,
2132                                     Stage.boundary_values,
2133                                     Xmom.boundary_values,
2134                                     Ymom.boundary_values,
2135                                     Stage.explicit_update,
2136                                     Xmom.explicit_update,
2137                                     Ymom.explicit_update,
2138                                     already_computed_flux)
2139
2140
2141    Post conditions:
2142      domain.explicit_update is reset to computed flux values
2143      domain.timestep is set to the largest step satisfying all volumes.
2144
2145
2146  */
2147
2148
2149  PyArrayObject *neighbours, *neighbour_edges,
2150    *normals, *edgelengths, *radii, *areas,
2151    *tri_full_flag,
2152    *stage_edge_values,
2153    *xmom_edge_values,
2154    *ymom_edge_values,
2155    *bed_edge_values,
2156    *stage_boundary_values,
2157    *xmom_boundary_values,
2158    *ymom_boundary_values,
2159    *stage_explicit_update,
2160    *xmom_explicit_update,
2161    *ymom_explicit_update,
2162    *already_computed_flux; // Tracks whether the flux across an edge has already been computed
2163
2164
2165  // Local variables
2166  double timestep, max_speed, epsilon, g, H0;
2167  double normal[2], ql[3], qr[3], zl, zr;
2168  double edgeflux[3]; //Work arrays for summing up fluxes
2169
2170  int number_of_elements, k, i, m, n;
2171  int ki, nm=0, ki2; //Index shorthands
2172  static long call=1;
2173
2174
2175  // Convert Python arguments to C
2176  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOO",
2177                        &timestep,
2178                        &epsilon,
2179                        &H0,
2180                        &g,
2181                        &neighbours,
2182                        &neighbour_edges,
2183                        &normals,
2184                        &edgelengths, &radii, &areas,
2185                        &tri_full_flag,
2186                        &stage_edge_values,
2187                        &xmom_edge_values,
2188                        &ymom_edge_values,
2189                        &bed_edge_values,
2190                        &stage_boundary_values,
2191                        &xmom_boundary_values,
2192                        &ymom_boundary_values,
2193                        &stage_explicit_update,
2194                        &xmom_explicit_update,
2195                        &ymom_explicit_update,
2196                        &already_computed_flux)) {
2197    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
2198    return NULL;
2199  }
2200  number_of_elements = stage_edge_values -> dimensions[0];
2201  call++;//a static local variable to which already_computed_flux is compared
2202  //set explicit_update to zero for all conserved_quantities.
2203  //This assumes compute_fluxes called before forcing terms
2204  for (k=0; k<number_of_elements; k++) {
2205    ((double *) stage_explicit_update -> data)[k]=0.0;
2206    ((double *) xmom_explicit_update -> data)[k]=0.0;
2207    ((double *) ymom_explicit_update -> data)[k]=0.0;
2208  }
2209  //Loop through neighbours and compute edge flux for each
2210  for (k=0; k<number_of_elements; k++) {
2211    for (i=0; i<3; i++) {
2212      ki = k*3+i;
2213      if (((long *) already_computed_flux->data)[ki]==call)//we've already computed the flux across this edge
2214        continue;
2215      ql[0] = ((double *) stage_edge_values -> data)[ki];
2216      ql[1] = ((double *) xmom_edge_values -> data)[ki];
2217      ql[2] = ((double *) ymom_edge_values -> data)[ki];
2218      zl =    ((double *) bed_edge_values -> data)[ki];
2219
2220      //Quantities at neighbour on nearest face
2221      n = ((long *) neighbours -> data)[ki];
2222      if (n < 0) {
2223        m = -n-1; //Convert negative flag to index
2224        qr[0] = ((double *) stage_boundary_values -> data)[m];
2225        qr[1] = ((double *) xmom_boundary_values -> data)[m];
2226        qr[2] = ((double *) ymom_boundary_values -> data)[m];
2227        zr = zl; //Extend bed elevation to boundary
2228      } else {
2229        m = ((long *) neighbour_edges -> data)[ki];
2230        nm = n*3+m;
2231        qr[0] = ((double *) stage_edge_values -> data)[nm];
2232        qr[1] = ((double *) xmom_edge_values -> data)[nm];
2233        qr[2] = ((double *) ymom_edge_values -> data)[nm];
2234        zr =    ((double *) bed_edge_values -> data)[nm];
2235      }
2236      // Outward pointing normal vector
2237      // normal = domain.normals[k, 2*i:2*i+2]
2238      ki2 = 2*ki; //k*6 + i*2
2239      normal[0] = ((double *) normals -> data)[ki2];
2240      normal[1] = ((double *) normals -> data)[ki2+1];
2241      //Edge flux computation
2242      flux_function_kinetic(ql, qr, zl, zr,
2243                    normal[0], normal[1],
2244                    epsilon, H0, g,
2245                    edgeflux, &max_speed);
2246      //update triangle k
2247      ((long *) already_computed_flux->data)[ki]=call;
2248      ((double *) stage_explicit_update -> data)[k] -= edgeflux[0]*((double *) edgelengths -> data)[ki];
2249      ((double *) xmom_explicit_update -> data)[k] -= edgeflux[1]*((double *) edgelengths -> data)[ki];
2250      ((double *) ymom_explicit_update -> data)[k] -= edgeflux[2]*((double *) edgelengths -> data)[ki];
2251      //update the neighbour n
2252      if (n>=0){
2253        ((long *) already_computed_flux->data)[nm]=call;
2254        ((double *) stage_explicit_update -> data)[n] += edgeflux[0]*((double *) edgelengths -> data)[nm];
2255        ((double *) xmom_explicit_update -> data)[n] += edgeflux[1]*((double *) edgelengths -> data)[nm];
2256        ((double *) ymom_explicit_update -> data)[n] += edgeflux[2]*((double *) edgelengths -> data)[nm];
2257      }
2258      ///for (j=0; j<3; j++) {
2259        ///flux[j] -= edgeflux[j]*((double *) edgelengths -> data)[ki];
2260        ///}
2261        //Update timestep
2262        //timestep = min(timestep, domain.radii[k]/max_speed)
2263        //FIXME: SR Add parameter for CFL condition
2264    if ( ((long *) tri_full_flag -> data)[k] == 1) {
2265            if (max_speed > epsilon) {
2266                timestep = min(timestep, ((double *) radii -> data)[k]/max_speed);
2267                //maxspeed in flux_function is calculated as max(|u+a|,|u-a|)
2268                if (n>=0)
2269                    timestep = min(timestep, ((double *) radii -> data)[n]/max_speed);
2270            }
2271    }
2272    } // end for i
2273    //Normalise by area and store for when all conserved
2274    //quantities get updated
2275    ((double *) stage_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
2276    ((double *) xmom_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
2277    ((double *) ymom_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
2278  } //end for k
2279  return Py_BuildValue("d", timestep);
2280}
2281
2282PyObject *protect(PyObject *self, PyObject *args) {
2283  //
2284  //    protect(minimum_allowed_height, maximum_allowed_speed, wc, zc, xmomc, ymomc)
2285
2286
2287  PyArrayObject
2288  *wc,            //Stage at centroids
2289  *zc,            //Elevation at centroids
2290  *xmomc,         //Momentums at centroids
2291  *ymomc;
2292
2293
2294  int N;
2295  double minimum_allowed_height, maximum_allowed_speed, epsilon;
2296
2297  // Convert Python arguments to C
2298  if (!PyArg_ParseTuple(args, "dddOOOO",
2299                        &minimum_allowed_height,
2300                        &maximum_allowed_speed,
2301                        &epsilon,
2302                        &wc, &zc, &xmomc, &ymomc)) {
2303    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: protect could not parse input arguments");
2304    return NULL;
2305  } 
2306
2307  N = wc -> dimensions[0];
2308
2309  _protect(N,
2310           minimum_allowed_height,
2311           maximum_allowed_speed,
2312           epsilon,
2313           (double*) wc -> data,
2314           (double*) zc -> data,
2315           (double*) xmomc -> data,
2316           (double*) ymomc -> data);
2317
2318  return Py_BuildValue("");
2319}
2320
2321
2322
2323PyObject *balance_deep_and_shallow(PyObject *self, PyObject *args) {
2324  // Compute linear combination between stage as computed by
2325  // gradient-limiters limiting using w, and stage computed by
2326  // gradient-limiters limiting using h (h-limiter).
2327  // The former takes precedence when heights are large compared to the
2328  // bed slope while the latter takes precedence when heights are
2329  // relatively small.  Anything in between is computed as a balanced
2330  // linear combination in order to avoid numerical disturbances which
2331  // would otherwise appear as a result of hard switching between
2332  // modes.
2333  //
2334  //    balance_deep_and_shallow(beta_h, wc, zc, wv, zv,
2335  //                             xmomc, ymomc, xmomv, ymomv)
2336
2337
2338  PyArrayObject
2339    *wc,            //Stage at centroids
2340    *zc,            //Elevation at centroids
2341    *wv,            //Stage at vertices
2342    *zv,            //Elevation at vertices
2343    *hvbar,         //h-Limited depths at vertices
2344    *xmomc,         //Momentums at centroids and vertices
2345    *ymomc,
2346    *xmomv,
2347    *ymomv;
2348 
2349  PyObject *domain, *Tmp;
2350   
2351  double alpha_balance = 2.0;
2352  double H0, beta_h;
2353
2354  int N, tight_slope_limiters, use_centroid_velocities; //, err;
2355
2356  // Convert Python arguments to C
2357  if (!PyArg_ParseTuple(args, "OdOOOOOOOOO",
2358                        &domain,
2359                        &beta_h,
2360                        &wc, &zc, 
2361                        &wv, &zv, &hvbar,
2362                        &xmomc, &ymomc, &xmomv, &ymomv)) {
2363    PyErr_SetString(PyExc_RuntimeError, 
2364                    "shallow_water_ext.c: balance_deep_and_shallow could not parse input arguments");
2365    return NULL;
2366  } 
2367         
2368         
2369  // FIXME (Ole): I tested this without GetAttrString and got time down
2370  // marginally from 4.0s to 3.8s. Consider passing everything in
2371  // through ParseTuple and profile.
2372 
2373  // Pull out parameters
2374  Tmp = PyObject_GetAttrString(domain, "alpha_balance");
2375  if (!Tmp) {
2376    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not obtain object alpha_balance from domain");
2377    return NULL;
2378  } 
2379  alpha_balance = PyFloat_AsDouble(Tmp);
2380  Py_DECREF(Tmp);
2381
2382 
2383  Tmp = PyObject_GetAttrString(domain, "H0");
2384  if (!Tmp) {
2385    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not obtain object H0 from domain");
2386    return NULL;
2387  } 
2388  H0 = PyFloat_AsDouble(Tmp);
2389  Py_DECREF(Tmp);
2390
2391 
2392  Tmp = PyObject_GetAttrString(domain, "tight_slope_limiters");
2393  if (!Tmp) {
2394    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not obtain object tight_slope_limiters from domain");
2395    return NULL;
2396  } 
2397  tight_slope_limiters = PyInt_AsLong(Tmp);
2398  Py_DECREF(Tmp);
2399 
2400 
2401  Tmp = PyObject_GetAttrString(domain, "use_centroid_velocities");
2402  if (!Tmp) {
2403    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not obtain object use_centroid_velocities from domain");
2404    return NULL;
2405  } 
2406  use_centroid_velocities = PyInt_AsLong(Tmp);
2407  Py_DECREF(Tmp);
2408 
2409 
2410 
2411  N = wc -> dimensions[0];
2412  _balance_deep_and_shallow(N,
2413                            beta_h,
2414                            (double*) wc -> data,
2415                            (double*) zc -> data,
2416                            (double*) wv -> data,
2417                            (double*) zv -> data,
2418                            (double*) hvbar -> data,
2419                            (double*) xmomc -> data,
2420                            (double*) ymomc -> data,
2421                            (double*) xmomv -> data,
2422                            (double*) ymomv -> data,
2423                            H0,
2424                            (int) tight_slope_limiters,
2425                            (int) use_centroid_velocities,                         
2426                            alpha_balance);
2427
2428
2429  return Py_BuildValue("");
2430}
2431
2432
2433
2434PyObject *h_limiter(PyObject *self, PyObject *args) {
2435
2436  PyObject *domain, *Tmp;
2437  PyArrayObject
2438    *hv, *hc, //Depth at vertices and centroids
2439    *hvbar,   //Limited depth at vertices (return values)
2440    *neighbours;
2441
2442  int k, i, n, N, k3;
2443  int dimensions[2];
2444  double beta_h; // Safety factor (see config.py)
2445  double *hmin, *hmax, hn;
2446
2447  // Convert Python arguments to C
2448  if (!PyArg_ParseTuple(args, "OOO", &domain, &hc, &hv)) {
2449    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: h_limiter could not parse input arguments");
2450    return NULL;
2451  } 
2452 
2453  neighbours = get_consecutive_array(domain, "neighbours");
2454
2455  //Get safety factor beta_h
2456  Tmp = PyObject_GetAttrString(domain, "beta_h");
2457  if (!Tmp) {
2458    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: h_limiter could not obtain object beta_h from domain");
2459    return NULL;
2460  } 
2461  beta_h = PyFloat_AsDouble(Tmp);
2462
2463  Py_DECREF(Tmp);
2464
2465  N = hc -> dimensions[0];
2466
2467  //Create hvbar
2468  dimensions[0] = N;
2469  dimensions[1] = 3;
2470  hvbar = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
2471
2472
2473  //Find min and max of this and neighbour's centroid values
2474  hmin = malloc(N * sizeof(double));
2475  hmax = malloc(N * sizeof(double));
2476  for (k=0; k<N; k++) {
2477    k3=k*3;
2478
2479    hmin[k] = ((double*) hc -> data)[k];
2480    hmax[k] = hmin[k];
2481
2482    for (i=0; i<3; i++) {
2483      n = ((long*) neighbours -> data)[k3+i];
2484
2485      // Initialise hvbar with values from hv
2486      ((double*) hvbar -> data)[k3+i] = ((double*) hv -> data)[k3+i];
2487
2488      if (n >= 0) {
2489        hn = ((double*) hc -> data)[n]; // Neighbour's centroid value
2490
2491        hmin[k] = min(hmin[k], hn);
2492        hmax[k] = max(hmax[k], hn);
2493      }
2494    }
2495  }
2496
2497  // Call underlying standard routine
2498  _limit_old(N, beta_h, (double*) hc -> data, (double*) hvbar -> data, hmin, hmax);
2499
2500  // // //Py_DECREF(domain); //FIXME: Necessary?
2501  free(hmin);
2502  free(hmax);
2503
2504  // Return result using PyArray to avoid memory leak
2505  return PyArray_Return(hvbar);
2506}
2507
2508PyObject *h_limiter_sw(PyObject *self, PyObject *args) {
2509  //a faster version of h_limiter above
2510  PyObject *domain, *Tmp;
2511  PyArrayObject
2512    *hv, *hc, //Depth at vertices and centroids
2513    *hvbar,   //Limited depth at vertices (return values)
2514    *neighbours;
2515
2516  int k, i, N, k3,k0,k1,k2;
2517  int dimensions[2];
2518  double beta_h; //Safety factor (see config.py)
2519  double hmin, hmax, dh[3];
2520// Convert Python arguments to C
2521  if (!PyArg_ParseTuple(args, "OOO", &domain, &hc, &hv)) {
2522    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: h_limiter_sw could not parse input arguments");
2523    return NULL;
2524  } 
2525  neighbours = get_consecutive_array(domain, "neighbours");
2526
2527  //Get safety factor beta_h
2528  Tmp = PyObject_GetAttrString(domain, "beta_h");
2529  if (!Tmp) {
2530    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: h_limiter_sw could not obtain object beta_h from domain");
2531    return NULL;
2532  } 
2533  beta_h = PyFloat_AsDouble(Tmp);
2534
2535  Py_DECREF(Tmp);
2536
2537  N = hc -> dimensions[0];
2538
2539  //Create hvbar
2540  dimensions[0] = N;
2541  dimensions[1] = 3;
2542  hvbar = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
2543  for (k=0;k<N;k++){
2544    k3=k*3;
2545    //get the ids of the neighbours
2546    k0 = ((long*) neighbours -> data)[k3];
2547    k1 = ((long*) neighbours -> data)[k3+1];
2548    k2 = ((long*) neighbours -> data)[k3+2];
2549    //set hvbar provisionally
2550    for (i=0;i<3;i++){
2551      ((double*) hvbar -> data)[k3+i] = ((double*) hv -> data)[k3+i];
2552      dh[i]=((double*) hvbar -> data)[k3+i]-((double*) hc -> data)[k];
2553    }
2554    hmin=((double*) hc -> data)[k];
2555    hmax=hmin;
2556    if (k0>=0){
2557      hmin=min(hmin,((double*) hc -> data)[k0]);
2558      hmax=max(hmax,((double*) hc -> data)[k0]);
2559    }
2560    if (k1>=0){
2561      hmin=min(hmin,((double*) hc -> data)[k1]);
2562      hmax=max(hmax,((double*) hc -> data)[k1]);
2563    }
2564    if (k2>=0){
2565      hmin=min(hmin,((double*) hc -> data)[k2]);
2566      hmax=max(hmax,((double*) hc -> data)[k2]);
2567    }
2568    hmin-=((double*) hc -> data)[k];
2569    hmax-=((double*) hc -> data)[k];
2570    limit_gradient(dh,hmin,hmax,beta_h);
2571    for (i=0;i<3;i++)
2572      ((double*) hvbar -> data)[k3+i] = ((double*) hc -> data)[k]+dh[i];
2573  }
2574  return PyArray_Return(hvbar);
2575}
2576
2577PyObject *assign_windfield_values(PyObject *self, PyObject *args) {
2578  //
2579  //      assign_windfield_values(xmom_update, ymom_update,
2580  //                              s_vec, phi_vec, self.const)
2581
2582
2583
2584  PyArrayObject   //(one element per triangle)
2585  *s_vec,         //Speeds
2586  *phi_vec,       //Bearings
2587  *xmom_update,   //Momentum updates
2588  *ymom_update;
2589
2590
2591  int N;
2592  double cw;
2593
2594  // Convert Python arguments to C
2595  if (!PyArg_ParseTuple(args, "OOOOd",
2596                        &xmom_update,
2597                        &ymom_update,
2598                        &s_vec, &phi_vec,
2599                        &cw)) {
2600    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: assign_windfield_values could not parse input arguments");
2601    return NULL;
2602  } 
2603                       
2604
2605  N = xmom_update -> dimensions[0];
2606
2607  _assign_wind_field_values(N,
2608           (double*) xmom_update -> data,
2609           (double*) ymom_update -> data,
2610           (double*) s_vec -> data,
2611           (double*) phi_vec -> data,
2612           cw);
2613
2614  return Py_BuildValue("");
2615}
2616
2617
2618
2619
2620//-------------------------------
2621// Method table for python module
2622//-------------------------------
2623static struct PyMethodDef MethodTable[] = {
2624  /* The cast of the function is necessary since PyCFunction values
2625   * only take two PyObject* parameters, and rotate() takes
2626   * three.
2627   */
2628
2629  {"rotate", (PyCFunction)rotate, METH_VARARGS | METH_KEYWORDS, "Print out"},
2630  {"extrapolate_second_order_sw", extrapolate_second_order_sw, METH_VARARGS, "Print out"},
2631  {"compute_fluxes_ext_central", compute_fluxes_ext_central, METH_VARARGS, "Print out"},
2632  {"compute_fluxes_ext_central_new", compute_fluxes_ext_central_new, METH_VARARGS, "Print out"},
2633  {"compute_fluxes_ext_kinetic", compute_fluxes_ext_kinetic, METH_VARARGS, "Print out"},
2634  {"gravity", gravity, METH_VARARGS, "Print out"},
2635  {"manning_friction", manning_friction, METH_VARARGS, "Print out"},
2636  {"flux_function_central", flux_function_central, METH_VARARGS, "Print out"}, 
2637  {"balance_deep_and_shallow", balance_deep_and_shallow,
2638   METH_VARARGS, "Print out"},
2639  {"h_limiter", h_limiter,
2640   METH_VARARGS, "Print out"},
2641  {"h_limiter_sw", h_limiter_sw,
2642   METH_VARARGS, "Print out"},
2643  {"protect", protect, METH_VARARGS | METH_KEYWORDS, "Print out"},
2644  {"assign_windfield_values", assign_windfield_values,
2645   METH_VARARGS | METH_KEYWORDS, "Print out"},
2646  {NULL, NULL}
2647};
2648
2649// Module initialisation
2650void initshallow_water_ext(void){
2651  Py_InitModule("shallow_water_ext", MethodTable);
2652
2653  import_array(); // Necessary for handling of NumPY structures
2654}
Note: See TracBrowser for help on using the repository browser.