# source:branches/numpy/anuga/shallow_water/shallow_water_ext.c@6902

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

Back-merge from Numeric trunk to numpy branch.

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