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

Last change on this file since 5162 was 5162, checked in by steve, 16 years ago

Updated some methods for quantity. Looks like we can use old
limiting system with larger values of beta.

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