# source:anuga_core/source/anuga/shallow_water/shallow_water_ext.c@5297

Last change on this file since 5297 was 5297, checked in by ole, 15 years ago

Implemented option use_centroid_velocities as per ticket:262.
The run_profile script on a linux amd 64 (bogong) improved from 37s to
33s or about 10%.

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