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

Last change on this file since 4251 was 4251, checked in by ole, 17 years ago

Pulled speed computation out into a separate function that can be used by all flux functions.

File size: 63.0 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
23const double pi = 3.14159265358979;
24
25// Computational function for rotation
26int _rotate(double *q, double n1, double n2) {
27  /*Rotate the momentum component q (q[1], q[2])
28    from x,y coordinates to coordinates based on normal vector (n1, n2).
29
30    Result is returned in array 3x1 r
31    To rotate in opposite direction, call rotate with (q, n1, -n2)
32
33    Contents of q are changed by this function */
34
35
36  double q1, q2;
37
38  //Shorthands
39  q1 = q[1];  //uh momentum
40  q2 = q[2];  //vh momentum
41
42  //Rotate
43  q[1] =  n1*q1 + n2*q2;
44  q[2] = -n2*q1 + n1*q2;
45
46  return 0;
47}
48
49int find_qmin_and_qmax(double dq0, double dq1, double dq2, double *qmin, double *qmax){
50  //Considering the centroid of an FV triangle and the vertices of its auxiliary triangle, find
51  //qmin=min(q)-qc and qmax=max(q)-qc, where min(q) and max(q) are respectively min and max over the
52  //four values (at the centroid of the FV triangle and the auxiliary triangle vertices),
53  //and qc is the centroid
54  //dq0=q(vertex0)-q(centroid of FV triangle)
55  //dq1=q(vertex1)-q(vertex0)
56  //dq2=q(vertex2)-q(vertex0)
57  if (dq0>=0.0){
58    if (dq1>=dq2){
59      if (dq1>=0.0)
60        *qmax=dq0+dq1;
61      else
62        *qmax=dq0;
63      if ((*qmin=dq0+dq2)<0)
64        ;//qmin is already set to correct value
65      else
66        *qmin=0.0;
67    }
68    else{//dq1<dq2
69      if (dq2>0)
70        *qmax=dq0+dq2;
71      else
72        *qmax=dq0;
73      if ((*qmin=dq0+dq1)<0)
74        ;//qmin is the correct value
75      else
76        *qmin=0.0;
77    }
78  }
79  else{//dq0<0
80    if (dq1<=dq2){
81      if (dq1<0.0)
82        *qmin=dq0+dq1;
83      else
84        *qmin=dq0;
85      if ((*qmax=dq0+dq2)>0.0)
86        ;//qmax is already set to the correct value
87      else
88        *qmax=0.0;
89    }
90    else{//dq1>dq2
91      if (dq2<0.0)
92        *qmin=dq0+dq2;
93      else
94        *qmin=dq0;
95      if ((*qmax=dq0+dq1)>0.0)
96        ;//qmax is already set to the correct value
97      else
98        *qmax=0.0;
99    }
100  }
101  return 0;
102}
103
104int limit_gradient(double *dqv, double qmin, double qmax, double beta_w){
105  //given provisional jumps dqv from the FV triangle centroid to its vertices and
106  //jumps qmin (qmax) between the centroid of the FV triangle and the
107  //minimum (maximum) of the values at the centroid of the FV triangle and the auxiliary triangle vertices,
108  //calculate a multiplicative factor phi by which the provisional vertex jumps are to be limited
109  int i;
110  double r=1000.0, r0=1.0, phi=1.0;
111  static double TINY = 1.0e-100;//to avoid machine accuracy problems.
112  //Any provisional jump with magnitude < TINY does not contribute to the limiting process.
113  for (i=0;i<3;i++){
114    if (dqv[i]<-TINY)
115      r0=qmin/dqv[i];
116    if (dqv[i]>TINY)
117      r0=qmax/dqv[i];
118    r=min(r0,r);
119    //
120  }
121  phi=min(r*beta_w,1.0);
122  for (i=0;i<3;i++)
123    dqv[i]=dqv[i]*phi;
124  return 0;
125}
126
127
128// Function to obtain speed from momentum and depth.
129// This is used by flux functions
130// Input parameters uh and h may be modified by this function.
131double _compute_speed(double *uh,
132                      double *h,
133                      double epsilon,
134                      double h0) {
135
136  double u;
137
138  if (*h < epsilon) {
139    *h = 0.0;  //Could have been negative
140    u = 0.0;
141  } else {
142    u = *uh/(*h + h0/ *h);
143  }
144
145  return u;
146}
147
148// Computational function for flux computation (using stage w=z+h)
149int flux_function_central(double *q_left, double *q_right,
150                  double z_left, double z_right,
151                  double n1, double n2,
152                  double epsilon, double H0, double g,
153                  double *edgeflux, double *max_speed) {
154
155  /*Compute fluxes between volumes for the shallow water wave equation
156    cast in terms of the 'stage', w = h+z using
157    the 'central scheme' as described in
158
159    Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
160    Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
161    Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.
162
163    The implemented formula is given in equation (3.15) on page 714
164  */
165
166  int i;
167
168  double w_left, h_left, uh_left, vh_left, u_left;
169  double w_right, h_right, uh_right, vh_right, u_right;
170  double s_min, s_max, soundspeed_left, soundspeed_right;
171  double denom, z;
172  double q_left_copy[3], q_right_copy[3];
173  double flux_right[3], flux_left[3];
174
175  double h0 = H0*H0; //This ensures a good balance when h approaches H0.
176
177  //Copy conserved quantities to protect from modification
178  for (i=0; i<3; i++) {
179    q_left_copy[i] = q_left[i];
180    q_right_copy[i] = q_right[i];
181  }
182
183  //Align x- and y-momentum with x-axis
184  _rotate(q_left_copy, n1, n2);
185  _rotate(q_right_copy, n1, n2);
186
187  z = (z_left+z_right)/2; //Take average of field values
188
189  //Compute speeds in x-direction
190  w_left = q_left_copy[0];
191  h_left = w_left-z;
192  uh_left = q_left_copy[1];
193  u_left =_compute_speed(&uh_left, &h_left, epsilon, h0);
194
195  w_right = q_right_copy[0];
196  h_right = w_right-z;
197  uh_right = q_right_copy[1];
198  u_right =_compute_speed(&uh_right, &h_right, epsilon, h0);
199
200  //Momentum in y-direction
201  vh_left  = q_left_copy[2];
202  vh_right = q_right_copy[2];
203
204
205  //Maximal and minimal wave speeds
206  soundspeed_left  = sqrt(g*h_left);
207  soundspeed_right = sqrt(g*h_right);
208
209  s_max = max(u_left+soundspeed_left, u_right+soundspeed_right);
210  if (s_max < 0.0) s_max = 0.0;
211
212  s_min = min(u_left-soundspeed_left, u_right-soundspeed_right);
213  if (s_min > 0.0) s_min = 0.0;
214
215  //Flux formulas
216  flux_left[0] = u_left*h_left;
217  flux_left[1] = u_left*uh_left + 0.5*g*h_left*h_left;
218  flux_left[2] = u_left*vh_left;
219
220  flux_right[0] = u_right*h_right;
221  flux_right[1] = u_right*uh_right + 0.5*g*h_right*h_right;
222  flux_right[2] = u_right*vh_right;
223
224
225  //Flux computation
226  denom = s_max-s_min;
227  if (denom == 0.0) {
228    for (i=0; i<3; i++) edgeflux[i] = 0.0;
229    *max_speed = 0.0;
230  } else {
231    for (i=0; i<3; i++) {
232      edgeflux[i] = s_max*flux_left[i] - s_min*flux_right[i];
233      edgeflux[i] += s_max*s_min*(q_right_copy[i]-q_left_copy[i]);
234      edgeflux[i] /= denom;
235    }
236
237    //Maximal wavespeed
238    *max_speed = max(fabs(s_max), fabs(s_min));
239
240    //Rotate back
241    _rotate(edgeflux, n1, -n2);
242  }
243  return 0;
244}
245
246double erfcc(double x){
247    double z,t,result;
248
249    z=fabs(x);
250    t=1.0/(1.0+0.5*z);
251    result=t*exp(-z*z-1.26551223+t*(1.00002368+t*(.37409196+
252         t*(.09678418+t*(-.18628806+t*(.27886807+t*(-1.13520398+
253         t*(1.48851587+t*(-.82215223+t*.17087277)))))))));
254    if (x < 0.0) result = 2.0-result;
255
256    return result;
257    }
258
259
260
261// Computational function for flux computation (using stage w=z+h)
262// FIXME (Ole): Is this used anywhere??
263int flux_function_kinetic(double *q_left, double *q_right,
264                  double z_left, double z_right,
265                  double n1, double n2,
266                  double epsilon, double H0, double g,
267                  double *edgeflux, double *max_speed) {
268
269  /*Compute fluxes between volumes for the shallow water wave equation
270    cast in terms of the 'stage', w = h+z using
271    the 'central scheme' as described in
272
273    Zhang et. al., Advances in Water Resources, 26(6), 2003, 635-647.
274  */
275
276  int i;
277
278  double w_left, h_left, uh_left, vh_left, u_left, F_left;
279  double w_right, h_right, uh_right, vh_right, u_right, F_right;
280  double s_min, s_max, soundspeed_left, soundspeed_right;
281  double z;
282  double q_left_copy[3], q_right_copy[3];
283
284  double h0 = H0*H0; //This ensures a good balance when h approaches H0.
285
286  //Copy conserved quantities to protect from modification
287  for (i=0; i<3; i++) {
288    q_left_copy[i] = q_left[i];
289    q_right_copy[i] = q_right[i];
290  }
291
292  //Align x- and y-momentum with x-axis
293  _rotate(q_left_copy, n1, n2);
294  _rotate(q_right_copy, n1, n2);
295
296  z = (z_left+z_right)/2; //Take average of field values
297
298  //Compute speeds in x-direction
299  w_left = q_left_copy[0];
300  h_left = w_left-z;
301  uh_left = q_left_copy[1];
302  u_left =_compute_speed(&uh_left, &h_left, epsilon, h0);
303
304  w_right = q_right_copy[0];
305  h_right = w_right-z;
306  uh_right = q_right_copy[1];
307  u_right =_compute_speed(&uh_right, &h_right, epsilon, h0);
308
309
310  //Momentum in y-direction
311  vh_left  = q_left_copy[2];
312  vh_right = q_right_copy[2];
313
314
315  //Maximal and minimal wave speeds
316  soundspeed_left  = sqrt(g*h_left);
317  soundspeed_right = sqrt(g*h_right);
318
319  s_max = max(u_left+soundspeed_left, u_right+soundspeed_right);
320  if (s_max < 0.0) s_max = 0.0;
321
322  s_min = min(u_left-soundspeed_left, u_right-soundspeed_right);
323  if (s_min > 0.0) s_min = 0.0;
324
325
326  F_left  = 0.0;
327  F_right = 0.0;
328  if (h_left > 0.0) F_left = u_left/sqrt(g*h_left);
329  if (h_right > 0.0) F_right = u_right/sqrt(g*h_right);
330
331  for (i=0; i<3; i++) edgeflux[i] = 0.0;
332  *max_speed = 0.0;
333
334  edgeflux[0] = h_left*u_left/2.0*erfcc(-F_left) +  \
335          h_left*sqrt(g*h_left)/2.0/sqrt(pi)*exp(-(F_left*F_left)) + \
336          h_right*u_right/2.0*erfcc(F_right) -  \
337          h_right*sqrt(g*h_right)/2.0/sqrt(pi)*exp(-(F_right*F_right));
338
339  edgeflux[1] = (h_left*u_left*u_left + g/2.0*h_left*h_left)/2.0*erfcc(-F_left) + \
340          u_left*h_left*sqrt(g*h_left)/2.0/sqrt(pi)*exp(-(F_left*F_left)) + \
341          (h_right*u_right*u_right + g/2.0*h_right*h_right)/2.0*erfcc(F_right) -  \
342          u_right*h_right*sqrt(g*h_right)/2.0/sqrt(pi)*exp(-(F_right*F_right));
343
344  edgeflux[2] = vh_left*u_left/2.0*erfcc(-F_left) + \
345          vh_left*sqrt(g*h_left)/2.0/sqrt(pi)*exp(-(F_left*F_left)) + \
346          vh_right*u_right/2.0*erfcc(F_right) - \
347          vh_right*sqrt(g*h_right)/2.0/sqrt(pi)*exp(-(F_right*F_right));
348
349  //Maximal wavespeed
350  *max_speed = max(fabs(s_max), fabs(s_min));
351
352  //Rotate back
353  _rotate(edgeflux, n1, -n2);
354
355  return 0;
356}
357
358
359
360
361void _manning_friction(double g, double eps, int N,
362                       double* w, double* z,
363                       double* uh, double* vh,
364                       double* eta, double* xmom, double* ymom) {
365
366  int k;
367  double S, h;
368
369  for (k=0; k<N; k++) {
370    if (eta[k] > eps) {
371      h = w[k]-z[k];
372      if (h >= eps) {
373        S = -g * eta[k]*eta[k] * sqrt((uh[k]*uh[k] + vh[k]*vh[k]));
374        S /= pow(h, 7.0/3);      //Expensive (on Ole's home computer)
375        //S /= exp(7.0/3.0*log(h));      //seems to save about 15% over manning_friction
376        //S /= h*h*(1 + h/3.0 - h*h/9.0); //FIXME: Could use a Taylor expansion
377
378
379        //Update momentum
380        xmom[k] += S*uh[k];
381        ymom[k] += S*vh[k];
382      }
383    }
384  }
385}
386
387
388/*
389void _manning_friction_explicit(double g, double eps, int N,
390                       double* w, double* z,
391                       double* uh, double* vh,
392                       double* eta, double* xmom, double* ymom) {
393
394  int k;
395  double S, h;
396
397  for (k=0; k<N; k++) {
398    if (eta[k] > eps) {
399      h = w[k]-z[k];
400      if (h >= eps) {
401        S = -g * eta[k]*eta[k] * sqrt((uh[k]*uh[k] + vh[k]*vh[k]));
402        S /= pow(h, 7.0/3);      //Expensive (on Ole's home computer)
403        //S /= exp(7.0/3.0*log(h));      //seems to save about 15% over manning_friction
404        //S /= h*h*(1 + h/3.0 - h*h/9.0); //FIXME: Could use a Taylor expansion
405
406
407        //Update momentum
408        xmom[k] += S*uh[k];
409        ymom[k] += S*vh[k];
410      }
411    }
412  }
413}
414*/
415
416int _balance_deep_and_shallow(int N,
417                              double* wc,
418                              double* zc,
419                              double* hc,
420                              double* wv,
421                              double* zv,
422                              double* hv,
423                              double* hvbar,
424                              double* xmomc,
425                              double* ymomc,
426                              double* xmomv,
427                              double* ymomv,
428                              double H0,
429                              int limit2007,
430                              double alpha_balance) {
431
432  int k, k3, i;
433  double dz, hmin, alpha, h_diff;
434
435  //Compute linear combination between w-limited stages and
436  //h-limited stages close to the bed elevation.
437
438  for (k=0; k<N; k++) {
439    // Compute maximal variation in bed elevation
440    // This quantitiy is
441    //     dz = max_i abs(z_i - z_c)
442    // and it is independent of dimension
443    // In the 1d case zc = (z0+z1)/2
444    // In the 2d case zc = (z0+z1+z2)/3
445
446    k3 = 3*k;
447
448    //FIXME: Try with this one precomputed
449    dz = 0.0;
450    hmin = hv[k3];
451    for (i=0; i<3; i++) {
452      if (limit2007 == 0) {
453        dz = max(dz, fabs(zv[k3+i]-zc[k]));
454      }
455
456      hmin = min(hmin, hv[k3+i]);
457    }
458
459
460    //Create alpha in [0,1], where alpha==0 means using the h-limited
461    //stage and alpha==1 means using the w-limited stage as
462    //computed by the gradient limiter (both 1st or 2nd order)
463
464
465    if (limit2007 == 0) {
466      //If hmin > dz/alpha_balance then alpha = 1 and the bed will have no
467      //effect
468      //If hmin < 0 then alpha = 0 reverting to constant height above bed.
469      //The parameter alpha_balance==2 by default
470
471
472      if (dz > 0.0) {
473        alpha = max( min( alpha_balance*hmin/dz, 1.0), 0.0 );
474      } else {
475        alpha = 1.0;  //Flat bed
476      }
477      //printf("Using old style limiter\n");
478
479    } else {
480
481      // 2007 Balanced Limiter
482
483      // Make alpha as large as possible but still ensure that
484      // final depth is positive
485
486      if (hmin < H0) {
487        alpha = 1.0;
488        for (i=0; i<3; i++) {
489
490          h_diff = hvbar[k3+i] - hv[k3+i];
491
492          if (h_diff <= 0) {
493            // Deep water triangle is further away from bed than
494            // shallow water (hbar < h). Any alpha will do
495
496          } else {
497            // Denominator is positive which means that we need some of the
498            // h-limited stage.
499
500            alpha = min(alpha, (hvbar[k3+i] - H0)/h_diff);
501          }
502        }
503
504        // Ensure alpha in [0,1]
505        if (alpha>1.0) alpha=1.0;
506        if (alpha<0.0) alpha=0.0;
507
508      } else {
509        // Use w-limited stage exclusively
510        alpha = 1.0;
511      }
512    }
513
514
515
516    //printf("k=%d, hmin=%.2f, dz=%.2f, alpha=%.2f, alpha_balance=%.2f\n",
517    //     k, hmin, dz, alpha, alpha_balance);
518
519    //printf("dz = %.3f, alpha = %.8f\n", dz, alpha);
520
521    //  Let
522    //
523    //    wvi be the w-limited stage (wvi = zvi + hvi)
524    //    wvi- be the h-limited state (wvi- = zvi + hvi-)
525    //
526    //
527    //  where i=0,1,2 denotes the vertex ids
528    //
529    //  Weighted balance between w-limited and h-limited stage is
530    //
531    //    wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi)
532    //
533    //  It follows that the updated wvi is
534    //    wvi := zvi + (1-alpha)*hvi- + alpha*hvi
535    //
536    //   Momentum is balanced between constant and limited
537
538    if (alpha < 1) {
539      for (i=0; i<3; i++) {
540        wv[k3+i] = zv[k3+i] + (1-alpha)*hvbar[k3+i] + alpha*hv[k3+i];
541
542        //Update momentum as a linear combination of
543        //xmomc and ymomc (shallow) and momentum
544        //from extrapolator xmomv and ymomv (deep).
545        //FIXME (Ole): Is this really needed?
546        xmomv[k3+i] = (1-alpha)*xmomc[k] + alpha*xmomv[k3+i];
547        ymomv[k3+i] = (1-alpha)*ymomc[k] + alpha*ymomv[k3+i];
548      }
549    }
550  }
551  return 0;
552}
553
554
555
556int _protect(int N,
557             double minimum_allowed_height,
558             double maximum_allowed_speed,
559             double epsilon,
560             double* wc,
561             double* zc,
562             double* xmomc,
563             double* ymomc) {
564
565  int k;
566  double hc;
567  double u, v, reduced_speed;
568
569  //Protect against initesimal and negative heights
570  for (k=0; k<N; k++) {
571    hc = wc[k] - zc[k];
572
573    if (hc < minimum_allowed_height) {
574
575      //Old code: Set momentum to zero and ensure h is non negative
576      //xmomc[k] = 0.0;
577      //ymomc[k] = 0.0;
578      //if (hc <= 0.0) wc[k] = zc[k];
579
580
581      //New code: Adjust momentum to guarantee speeds are physical
582      //          ensure h is non negative
583      //FIXME (Ole): This is only implemented in this C extension and
584      //             has no Python equivalent
585
586      if (hc <= 0.0) {
587        wc[k] = zc[k];
588        xmomc[k] = 0.0;
589        ymomc[k] = 0.0;
590      } else {
591        //Reduce excessive speeds derived from division by small hc
592        //FIXME (Ole): This may be unnecessary with new slope limiters
593        //in effect.
594
595        u = xmomc[k]/hc;
596        if (fabs(u) > maximum_allowed_speed) {
597          reduced_speed = maximum_allowed_speed * u/fabs(u);
598          //printf("Speed (u) has been reduced from %.3f to %.3f\n",
599          //     u, reduced_speed);
600          xmomc[k] = reduced_speed * hc;
601        }
602
603        v = ymomc[k]/hc;
604        if (fabs(v) > maximum_allowed_speed) {
605          reduced_speed = maximum_allowed_speed * v/fabs(v);
606          //printf("Speed (v) has been reduced from %.3f to %.3f\n",
607          //     v, reduced_speed);
608          ymomc[k] = reduced_speed * hc;
609        }
610      }
611    }
612  }
613  return 0;
614}
615
616
617
618
619int _assign_wind_field_values(int N,
620                              double* xmom_update,
621                              double* ymom_update,
622                              double* s_vec,
623                              double* phi_vec,
624                              double cw) {
625
626  //Assign windfield values to momentum updates
627
628  int k;
629  double S, s, phi, u, v;
630
631  for (k=0; k<N; k++) {
632
633    s = s_vec[k];
634    phi = phi_vec[k];
635
636    //Convert to radians
637    phi = phi*pi/180;
638
639    //Compute velocity vector (u, v)
640    u = s*cos(phi);
641    v = s*sin(phi);
642
643    //Compute wind stress
644    S = cw * sqrt(u*u + v*v);
645    xmom_update[k] += S*u;
646    ymom_update[k] += S*v;
647  }
648  return 0;
649}
650
651
652
653///////////////////////////////////////////////////////////////////
654// Gateways to Python
655
656PyObject *gravity(PyObject *self, PyObject *args) {
657  //
658  //  gravity(g, h, v, x, xmom, ymom)
659  //
660
661
662  PyArrayObject *h, *v, *x, *xmom, *ymom;
663  int k, i, N, k3, k6;
664  double g, avg_h, zx, zy;
665  double x0, y0, x1, y1, x2, y2, z0, z1, z2;
666
667  if (!PyArg_ParseTuple(args, "dOOOOO",
668                        &g, &h, &v, &x,
669                        &xmom, &ymom)) {
670    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: gravity could not parse input arguments");
671    return NULL;
672  }
673
674  N = h -> dimensions[0];
675  for (k=0; k<N; k++) {
676    k3 = 3*k;  // base index
677    k6 = 6*k;  // base index
678
679    avg_h = 0.0;
680    for (i=0; i<3; i++) {
681      avg_h += ((double *) h -> data)[k3+i];
682    }
683    avg_h /= 3;
684
685
686    //Compute bed slope
687    x0 = ((double*) x -> data)[k6 + 0];
688    y0 = ((double*) x -> data)[k6 + 1];
689    x1 = ((double*) x -> data)[k6 + 2];
690    y1 = ((double*) x -> data)[k6 + 3];
691    x2 = ((double*) x -> data)[k6 + 4];
692    y2 = ((double*) x -> data)[k6 + 5];
693
694
695    z0 = ((double*) v -> data)[k3 + 0];
696    z1 = ((double*) v -> data)[k3 + 1];
697    z2 = ((double*) v -> data)[k3 + 2];
698
699    _gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2, &zx, &zy);
700
701    //Update momentum
702    ((double*) xmom -> data)[k] += -g*zx*avg_h;
703    ((double*) ymom -> data)[k] += -g*zy*avg_h;
704  }
705
706  return Py_BuildValue("");
707}
708
709
710PyObject *manning_friction(PyObject *self, PyObject *args) {
711  //
712  // manning_friction(g, eps, h, uh, vh, eta, xmom_update, ymom_update)
713  //
714
715
716  PyArrayObject *w, *z, *uh, *vh, *eta, *xmom, *ymom;
717  int N;
718  double g, eps;
719
720  if (!PyArg_ParseTuple(args, "ddOOOOOOO",
721                        &g, &eps, &w, &z, &uh, &vh, &eta,
722                        &xmom, &ymom)) {
723    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: manning_friction could not parse input arguments");
724    return NULL;
725  }
726
727
728  N = w -> dimensions[0];
729  _manning_friction(g, eps, N,
730                    (double*) w -> data,
731                    (double*) z -> data,
732                    (double*) uh -> data,
733                    (double*) vh -> data,
734                    (double*) eta -> data,
735                    (double*) xmom -> data,
736                    (double*) ymom -> data);
737
738  return Py_BuildValue("");
739}
740
741
742/*
743PyObject *manning_friction_explicit(PyObject *self, PyObject *args) {
744  //
745  // manning_friction_explicit(g, eps, h, uh, vh, eta, xmom_update, ymom_update)
746  //
747
748
749  PyArrayObject *w, *z, *uh, *vh, *eta, *xmom, *ymom;
750  int N;
751  double g, eps;
752
753  if (!PyArg_ParseTuple(args, "ddOOOOOOO",
754                        &g, &eps, &w, &z, &uh, &vh, &eta,
755                        &xmom, &ymom))
756    return NULL;
757
758  N = w -> dimensions[0];
759  _manning_friction_explicit(g, eps, N,
760                    (double*) w -> data,
761                    (double*) z -> data,
762                    (double*) uh -> data,
763                    (double*) vh -> data,
764                    (double*) eta -> data,
765                    (double*) xmom -> data,
766                    (double*) ymom -> data);
767
768  return Py_BuildValue("");
769}
770*/
771
772PyObject *extrapolate_second_order_sw(PyObject *self, PyObject *args) {
773  /*Compute the vertex values based on a linear reconstruction on each triangle
774    These values are calculated as follows:
775    1) For each triangle not adjacent to a boundary, we consider the auxiliary triangle
776    formed by the centroids of its three neighbours.
777    2) For each conserved quantity, we integrate around the auxiliary triangle's boundary the product
778    of the quantity and the outward normal vector. Dividing by the triangle area gives (a,b), the average
779    of the vector (q_x,q_y) on the auxiliary triangle. We suppose that the linear reconstruction on the
780    original triangle has gradient (a,b).
781    3) Provisional vertex junmps dqv[0,1,2] are computed and these are then limited by calling the functions
782    find_qmin_and_qmax and limit_gradient
783
784    Python call:
785    extrapolate_second_order_sw(domain.surrogate_neighbours,
786                                domain.number_of_boundaries
787                                domain.centroid_coordinates,
788                                Stage.centroid_values
789                                Xmom.centroid_values
790                                Ymom.centroid_values
791                                domain.vertex_coordinates,
792                                Stage.vertex_values,
793                                Xmom.vertex_values,
794                                Ymom.vertex_values)
795
796    Post conditions:
797            The vertices of each triangle have values from a limited linear reconstruction
798            based on centroid values
799
800  */
801  PyArrayObject *surrogate_neighbours,
802    *number_of_boundaries,
803    *centroid_coordinates,
804    *stage_centroid_values,
805    *xmom_centroid_values,
806    *ymom_centroid_values,
807        *elevation_centroid_values,
808    *vertex_coordinates,
809    *stage_vertex_values,
810    *xmom_vertex_values,
811    *ymom_vertex_values,
812        *elevation_vertex_values;
813  PyObject *domain, *Tmp;
814  //Local variables
815  double a, b;//gradient vector, not stored but used to calculate vertex values from centroids
816  int number_of_elements,k,k0,k1,k2,k3,k6,coord_index,i;
817  double x,y,x0,y0,x1,y1,x2,y2,xv0,yv0,xv1,yv1,xv2,yv2;//vertices of the auxiliary triangle
818  double dx1,dx2,dy1,dy2,dxv0,dxv1,dxv2,dyv0,dyv1,dyv2,dq0,dq1,dq2,area2;
819  double dqv[3], qmin, qmax, hmin;
820  double hc, h0, h1, h2;
821  double beta_w, beta_w_dry, beta_uh, beta_uh_dry, beta_vh, beta_vh_dry, beta_tmp;
822  double minimum_allowed_height;
823  //provisional jumps from centroids to v'tices and safety factor re limiting
824  //by which these jumps are limited
825  // Convert Python arguments to C
826  if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOO",
827                        &domain,
828                        &surrogate_neighbours,
829                        &number_of_boundaries,
830                        &centroid_coordinates,
831                        &stage_centroid_values,
832                        &xmom_centroid_values,
833                        &ymom_centroid_values,
834                        &elevation_centroid_values,
835                        &vertex_coordinates,
836                        &stage_vertex_values,
837                        &xmom_vertex_values,
838                        &ymom_vertex_values,
839                        &elevation_vertex_values)) {
840    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
841    return NULL;
842  }
843
844  //get the safety factor beta_w, set in the config.py file. This is used in the limiting process
845  Tmp = PyObject_GetAttrString(domain, "beta_w");
846  if (!Tmp) {
847    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_w from domain");
848    return NULL;
849  }
850  beta_w = PyFloat_AsDouble(Tmp);
851  Py_DECREF(Tmp);
852
853  Tmp = PyObject_GetAttrString(domain, "beta_w_dry");
854  if (!Tmp) {
855    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_w_dry from domain");
856    return NULL;
857  }
858  beta_w_dry = PyFloat_AsDouble(Tmp);
859  Py_DECREF(Tmp);
860
861  Tmp = PyObject_GetAttrString(domain, "beta_uh");
862  if (!Tmp) {
863    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_uh from domain");
864    return NULL;
865  }
866  beta_uh = PyFloat_AsDouble(Tmp);
867  Py_DECREF(Tmp);
868
869  Tmp = PyObject_GetAttrString(domain, "beta_uh_dry");
870  if (!Tmp) {
871    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_uh_dry from domain");
872    return NULL;
873  }
874  beta_uh_dry = PyFloat_AsDouble(Tmp);
875  Py_DECREF(Tmp);
876
877  Tmp = PyObject_GetAttrString(domain, "beta_vh");
878  if (!Tmp) {
879    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_vh from domain");
880    return NULL;
881  }
882  beta_vh = PyFloat_AsDouble(Tmp);
883  Py_DECREF(Tmp);
884
885  Tmp = PyObject_GetAttrString(domain, "beta_vh_dry");
886  if (!Tmp) {
887    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_vh_dry from domain");
888    return NULL;
889  }
890  beta_vh_dry = PyFloat_AsDouble(Tmp);
891  Py_DECREF(Tmp);
892
893  Tmp = PyObject_GetAttrString(domain, "minimum_allowed_height");
894  if (!Tmp) {
895    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object minimum_allowed_heigt");
896    return NULL;
897  }
898  minimum_allowed_height = PyFloat_AsDouble(Tmp);
899  Py_DECREF(Tmp);
900
901  number_of_elements = stage_centroid_values -> dimensions[0];
902  for (k=0; k<number_of_elements; k++) {
903    k3=k*3;
904    k6=k*6;
905
906    if (((long *) number_of_boundaries->data)[k]==3){/*no neighbours, set gradient on the triangle to zero*/
907      ((double *) stage_vertex_values->data)[k3]=((double *)stage_centroid_values->data)[k];
908      ((double *) stage_vertex_values->data)[k3+1]=((double *)stage_centroid_values->data)[k];
909      ((double *) stage_vertex_values->data)[k3+2]=((double *)stage_centroid_values->data)[k];
910      ((double *) xmom_vertex_values->data)[k3]=((double *)xmom_centroid_values->data)[k];
911      ((double *) xmom_vertex_values->data)[k3+1]=((double *)xmom_centroid_values->data)[k];
912      ((double *) xmom_vertex_values->data)[k3+2]=((double *)xmom_centroid_values->data)[k];
913      ((double *) ymom_vertex_values->data)[k3]=((double *)ymom_centroid_values->data)[k];
914      ((double *) ymom_vertex_values->data)[k3+1]=((double *)ymom_centroid_values->data)[k];
915      ((double *) ymom_vertex_values->data)[k3+2]=((double *)ymom_centroid_values->data)[k];
916      continue;
917    }
918    else {//we will need centroid coordinates and vertex coordinates of the triangle
919      //get the vertex coordinates of the FV triangle
920      xv0=((double *)vertex_coordinates->data)[k6]; yv0=((double *)vertex_coordinates->data)[k6+1];
921      xv1=((double *)vertex_coordinates->data)[k6+2]; yv1=((double *)vertex_coordinates->data)[k6+3];
922      xv2=((double *)vertex_coordinates->data)[k6+4]; yv2=((double *)vertex_coordinates->data)[k6+5];
923      //get the centroid coordinates of the FV triangle
924      coord_index=2*k;
925      x=((double *)centroid_coordinates->data)[coord_index];
926      y=((double *)centroid_coordinates->data)[coord_index+1];
927      //store x- and y- differentials for the vertices of the FV triangle relative to the centroid
928      dxv0=xv0-x; dxv1=xv1-x; dxv2=xv2-x;
929      dyv0=yv0-y; dyv1=yv1-y; dyv2=yv2-y;
930    }
931    if (((long *)number_of_boundaries->data)[k]<=1){
932      //if no boundaries, auxiliary triangle is formed from the centroids of the three neighbours
933      //if one boundary, auxiliary triangle is formed from this centroid and its two neighbours
934      k0=((long *)surrogate_neighbours->data)[k3];
935      k1=((long *)surrogate_neighbours->data)[k3+1];
936      k2=((long *)surrogate_neighbours->data)[k3+2];
937      //get the auxiliary triangle's vertex coordinates (really the centroids of neighbouring triangles)
938      coord_index=2*k0;
939      x0=((double *)centroid_coordinates->data)[coord_index];
940      y0=((double *)centroid_coordinates->data)[coord_index+1];
941      coord_index=2*k1;
942      x1=((double *)centroid_coordinates->data)[coord_index];
943      y1=((double *)centroid_coordinates->data)[coord_index+1];
944      coord_index=2*k2;
945      x2=((double *)centroid_coordinates->data)[coord_index];
946      y2=((double *)centroid_coordinates->data)[coord_index+1];
947      //store x- and y- differentials for the vertices of the auxiliary triangle
948      dx1=x1-x0; dx2=x2-x0;
949      dy1=y1-y0; dy2=y2-y0;
950      //calculate 2*area of the auxiliary triangle
951      area2 = dy2*dx1 - dy1*dx2;//the triangle is guaranteed to be counter-clockwise
952      //If the mesh is 'weird' near the boundary, the trianlge might be flat or clockwise:
953      if (area2<=0) {
954        PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: negative triangle area encountered");
955        return NULL;
956      }
957
958
959      //### Calculate heights of neighbouring cells
960      hc = ((double *)stage_centroid_values->data)[k]  - ((double *)elevation_centroid_values->data)[k];
961      h0 = ((double *)stage_centroid_values->data)[k0] - ((double *)elevation_centroid_values->data)[k0];
962      h1 = ((double *)stage_centroid_values->data)[k1] - ((double *)elevation_centroid_values->data)[k1];
963      h2 = ((double *)stage_centroid_values->data)[k2] - ((double *)elevation_centroid_values->data)[k2];
964      hmin = min(hc,min(h0,min(h1,h2)));
965
966      //### stage ###
967      //calculate the difference between vertex 0 of the auxiliary triangle and the FV triangle centroid
968      dq0=((double *)stage_centroid_values->data)[k0]-((double *)stage_centroid_values->data)[k];
969      //calculate differentials between the vertices of the auxiliary triangle
970      dq1=((double *)stage_centroid_values->data)[k1]-((double *)stage_centroid_values->data)[k0];
971      dq2=((double *)stage_centroid_values->data)[k2]-((double *)stage_centroid_values->data)[k0];
972      //calculate the gradient of stage on the auxiliary triangle
973      a = dy2*dq1 - dy1*dq2;
974      a /= area2;
975      b = dx1*dq2 - dx2*dq1;
976      b /= area2;
977      //calculate provisional jumps in stage from the centroid of the FV tri to its vertices, to be limited
978      dqv[0]=a*dxv0+b*dyv0;
979      dqv[1]=a*dxv1+b*dyv1;
980      dqv[2]=a*dxv2+b*dyv2;
981      //now we want to find min and max of the centroid and the vertices of the auxiliary triangle
982      //and compute jumps from the centroid to the min and max
983      find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
984      // Playing with dry wet interface
985      hmin = qmin;
986      beta_tmp = beta_w;
987      if (hmin<minimum_allowed_height)
988        beta_tmp = beta_w_dry;
990      for (i=0;i<3;i++)
991        ((double *)stage_vertex_values->data)[k3+i]=((double *)stage_centroid_values->data)[k]+dqv[i];
992
993      //### xmom ###
994      //calculate the difference between vertex 0 of the auxiliary triangle and the FV triangle centroid
995      dq0=((double *)xmom_centroid_values->data)[k0]-((double *)xmom_centroid_values->data)[k];
996      //calculate differentials between the vertices of the auxiliary triangle
997      dq1=((double *)xmom_centroid_values->data)[k1]-((double *)xmom_centroid_values->data)[k0];
998      dq2=((double *)xmom_centroid_values->data)[k2]-((double *)xmom_centroid_values->data)[k0];
999      //calculate the gradient of xmom on the auxiliary triangle
1000      a = dy2*dq1 - dy1*dq2;
1001      a /= area2;
1002      b = dx1*dq2 - dx2*dq1;
1003      b /= area2;
1004      //calculate provisional jumps in stage from the centroid of the FV tri to its vertices, to be limited
1005      dqv[0]=a*dxv0+b*dyv0;
1006      dqv[1]=a*dxv1+b*dyv1;
1007      dqv[2]=a*dxv2+b*dyv2;
1008      //now we want to find min and max of the centroid and the vertices of the auxiliary triangle
1009      //and compute jumps from the centroid to the min and max
1010      find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
1011      beta_tmp = beta_uh;
1012      if (hmin<minimum_allowed_height)
1013        beta_tmp = beta_uh_dry;
1015      for (i=0;i<3;i++)
1016        ((double *)xmom_vertex_values->data)[k3+i]=((double *)xmom_centroid_values->data)[k]+dqv[i];
1017
1018      //### ymom ###
1019      //calculate the difference between vertex 0 of the auxiliary triangle and the FV triangle centroid
1020      dq0=((double *)ymom_centroid_values->data)[k0]-((double *)ymom_centroid_values->data)[k];
1021      //calculate differentials between the vertices of the auxiliary triangle
1022      dq1=((double *)ymom_centroid_values->data)[k1]-((double *)ymom_centroid_values->data)[k0];
1023      dq2=((double *)ymom_centroid_values->data)[k2]-((double *)ymom_centroid_values->data)[k0];
1024      //calculate the gradient of xmom on the auxiliary triangle
1025      a = dy2*dq1 - dy1*dq2;
1026      a /= area2;
1027      b = dx1*dq2 - dx2*dq1;
1028      b /= area2;
1029      //calculate provisional jumps in stage from the centroid of the FV tri to its vertices, to be limited
1030      dqv[0]=a*dxv0+b*dyv0;
1031      dqv[1]=a*dxv1+b*dyv1;
1032      dqv[2]=a*dxv2+b*dyv2;
1033      //now we want to find min and max of the centroid and the vertices of the auxiliary triangle
1034      //and compute jumps from the centroid to the min and max
1035      find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
1036      beta_tmp = beta_vh;
1037      if (hmin<minimum_allowed_height)
1038        beta_tmp = beta_vh_dry;
1040      for (i=0;i<3;i++)
1041        ((double *)ymom_vertex_values->data)[k3+i]=((double *)ymom_centroid_values->data)[k]+dqv[i];
1042    }//if (number_of_boundaries[k]<=1)
1043    else{//number_of_boundaries==2
1044      //one internal neighbour and gradient is in direction of the neighbour's centroid
1045      //find the only internal neighbour
1046      for (k2=k3;k2<k3+3;k2++){//k2 just indexes the edges of triangle k
1047        if (((long *)surrogate_neighbours->data)[k2]!=k)//find internal neighbour of triabngle k
1048          break;
1049      }
1050      if ((k2==k3+3)) {//if we didn't find an internal neighbour
1051        PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: Internal neighbour not found");
1052        return NULL;//error
1053      }
1054
1055      k1=((long *)surrogate_neighbours->data)[k2];
1056      //the coordinates of the triangle are already (x,y). Get centroid of the neighbour (x1,y1)
1057      coord_index=2*k1;
1058      x1=((double *)centroid_coordinates->data)[coord_index];
1059      y1=((double *)centroid_coordinates->data)[coord_index+1];
1060      //compute x- and y- distances between the centroid of the FV triangle and that of its neighbour
1061      dx1=x1-x; dy1=y1-y;
1062      //set area2 as the square of the distance
1063      area2=dx1*dx1+dy1*dy1;
1064      //set dx2=(x1-x0)/((x1-x0)^2+(y1-y0)^2) and dy2=(y1-y0)/((x1-x0)^2+(y1-y0)^2) which
1065      //respectively correspond to the x- and y- gradients of the conserved quantities
1066      dx2=1.0/area2;
1067      dy2=dx2*dy1;
1068      dx2*=dx1;
1069
1070      //## stage ###
1071      //compute differentials
1072      dq1=((double *)stage_centroid_values->data)[k1]-((double *)stage_centroid_values->data)[k];
1073      //calculate the gradient between the centroid of the FV triangle and that of its neighbour
1074      a=dq1*dx2;
1075      b=dq1*dy2;
1076      //calculate provisional vertex jumps, to be limited
1077      dqv[0]=a*dxv0+b*dyv0;
1078      dqv[1]=a*dxv1+b*dyv1;
1079      dqv[2]=a*dxv2+b*dyv2;
1080      //now limit the jumps
1081      if (dq1>=0.0){
1082        qmin=0.0;
1083        qmax=dq1;
1084      }
1085      else{
1086        qmin=dq1;
1087        qmax=0.0;
1088      }
1089
1090
1092      for (i=0;i<3;i++)
1093        ((double *)stage_vertex_values->data)[k3+i]=((double *)stage_centroid_values->data)[k]+dqv[i];
1094
1095      //## xmom ###
1096      //compute differentials
1097      dq1=((double *)xmom_centroid_values->data)[k1]-((double *)xmom_centroid_values->data)[k];
1098      //calculate the gradient between the centroid of the FV triangle and that of its neighbour
1099      a=dq1*dx2;
1100      b=dq1*dy2;
1101      //calculate provisional vertex jumps, to be limited
1102      dqv[0]=a*dxv0+b*dyv0;
1103      dqv[1]=a*dxv1+b*dyv1;
1104      dqv[2]=a*dxv2+b*dyv2;
1105      //now limit the jumps
1106      if (dq1>=0.0){
1107        qmin=0.0;
1108        qmax=dq1;
1109      }
1110      else{
1111        qmin=dq1;
1112        qmax=0.0;
1113      }
1115      for (i=0;i<3;i++)
1116        ((double *)xmom_vertex_values->data)[k3+i]=((double *)xmom_centroid_values->data)[k]+dqv[i];
1117
1118      //## ymom ###
1119      //compute differentials
1120      dq1=((double *)ymom_centroid_values->data)[k1]-((double *)ymom_centroid_values->data)[k];
1121      //calculate the gradient between the centroid of the FV triangle and that of its neighbour
1122      a=dq1*dx2;
1123      b=dq1*dy2;
1124      //calculate provisional vertex jumps, to be limited
1125      dqv[0]=a*dxv0+b*dyv0;
1126      dqv[1]=a*dxv1+b*dyv1;
1127      dqv[2]=a*dxv2+b*dyv2;
1128      //now limit the jumps
1129      if (dq1>=0.0){
1130        qmin=0.0;
1131        qmax=dq1;
1132      }
1133      else{
1134        qmin=dq1;
1135        qmax=0.0;
1136      }
1138      for (i=0;i<3;i++)
1139        ((double *)ymom_vertex_values->data)[k3+i]=((double *)ymom_centroid_values->data)[k]+dqv[i];
1140    }//else [number_of_boudaries==2]
1141  }//for k=0 to number_of_elements-1
1142  return Py_BuildValue("");
1143}//extrapolate_second-order_sw
1144
1145
1146PyObject *rotate(PyObject *self, PyObject *args, PyObject *kwargs) {
1147  //
1148  // r = rotate(q, normal, direction=1)
1149  //
1150  // Where q is assumed to be a Float numeric array of length 3 and
1151  // normal a Float numeric array of length 2.
1152
1153
1154  PyObject *Q, *Normal;
1155  PyArrayObject *q, *r, *normal;
1156
1157  static char *argnames[] = {"q", "normal", "direction", NULL};
1158  int dimensions[1], i, direction=1;
1159  double n1, n2;
1160
1161  // Convert Python arguments to C
1162  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OO|i", argnames,
1163                                   &Q, &Normal, &direction)) {
1164    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: rotate could not parse input arguments");
1165    return NULL;
1166  }
1167
1168  //Input checks (convert sequences into numeric arrays)
1169  q = (PyArrayObject *)
1170    PyArray_ContiguousFromObject(Q, PyArray_DOUBLE, 0, 0);
1171  normal = (PyArrayObject *)
1172    PyArray_ContiguousFromObject(Normal, PyArray_DOUBLE, 0, 0);
1173
1174
1175  if (normal -> dimensions[0] != 2) {
1176    PyErr_SetString(PyExc_RuntimeError, "Normal vector must have 2 components");
1177    return NULL;
1178  }
1179
1180  //Allocate space for return vector r (don't DECREF)
1181  dimensions[0] = 3;
1182  r = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
1183
1184  //Copy
1185  for (i=0; i<3; i++) {
1186    ((double *) (r -> data))[i] = ((double *) (q -> data))[i];
1187  }
1188
1189  //Get normal and direction
1190  n1 = ((double *) normal -> data)[0];
1191  n2 = ((double *) normal -> data)[1];
1192  if (direction == -1) n2 = -n2;
1193
1194  //Rotate
1195  _rotate((double *) r -> data, n1, n2);
1196
1197  //Release numeric arrays
1198  Py_DECREF(q);
1199  Py_DECREF(normal);
1200
1201  //return result using PyArray to avoid memory leak
1202  return PyArray_Return(r);
1203}
1204
1205PyObject *compute_fluxes_ext_central(PyObject *self, PyObject *args) {
1206  /*Compute all fluxes and the timestep suitable for all volumes
1207    in domain.
1208
1209    Compute total flux for each conserved quantity using "flux_function_central"
1210
1211    Fluxes across each edge are scaled by edgelengths and summed up
1212    Resulting flux is then scaled by area and stored in
1213    explicit_update for each of the three conserved quantities
1214    stage, xmomentum and ymomentum
1215
1216    The maximal allowable speed computed by the flux_function for each volume
1217    is converted to a timestep that must not be exceeded. The minimum of
1218    those is computed as the next overall timestep.
1219
1220    Python call:
1221    domain.timestep = compute_fluxes(timestep,
1222                                     domain.epsilon,
1223                                     domain.H0,
1224                                     domain.g,
1225                                     domain.neighbours,
1226                                     domain.neighbour_edges,
1227                                     domain.normals,
1228                                     domain.edgelengths,
1230                                     domain.areas,
1231                                     tri_full_flag,
1232                                     Stage.edge_values,
1233                                     Xmom.edge_values,
1234                                     Ymom.edge_values,
1235                                     Bed.edge_values,
1236                                     Stage.boundary_values,
1237                                     Xmom.boundary_values,
1238                                     Ymom.boundary_values,
1239                                     Stage.explicit_update,
1240                                     Xmom.explicit_update,
1241                                     Ymom.explicit_update,
1243
1244
1245    Post conditions:
1246      domain.explicit_update is reset to computed flux values
1247      domain.timestep is set to the largest step satisfying all volumes.
1248
1249
1250  */
1251
1252
1253  PyArrayObject *neighbours, *neighbour_edges,
1254    *normals, *edgelengths, *radii, *areas,
1255    *tri_full_flag,
1256    *stage_edge_values,
1257    *xmom_edge_values,
1258    *ymom_edge_values,
1259    *bed_edge_values,
1260    *stage_boundary_values,
1261    *xmom_boundary_values,
1262    *ymom_boundary_values,
1263    *stage_explicit_update,
1264    *xmom_explicit_update,
1265    *ymom_explicit_update,
1266    *already_computed_flux,//tracks whether the flux across an edge has already been computed
1267    *max_speed_array; //Keeps track of max speeds for each triangle
1268
1269  //Local variables
1270  double timestep, max_speed, epsilon, g, H0;
1271  double normal[2], ql[3], qr[3], zl, zr;
1272  double edgeflux[3]; //Work arrays for summing up fluxes
1273
1274  int number_of_elements, k, i, m, n;
1275  int ki, nm=0, ki2; //Index shorthands
1276  static long call=1;
1277
1278
1279  // Convert Python arguments to C
1280  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOO",
1281                        &timestep,
1282                        &epsilon,
1283                        &H0,
1284                        &g,
1285                        &neighbours,
1286                        &neighbour_edges,
1287                        &normals,
1288                        &edgelengths, &radii, &areas,
1289                        &tri_full_flag,
1290                        &stage_edge_values,
1291                        &xmom_edge_values,
1292                        &ymom_edge_values,
1293                        &bed_edge_values,
1294                        &stage_boundary_values,
1295                        &xmom_boundary_values,
1296                        &ymom_boundary_values,
1297                        &stage_explicit_update,
1298                        &xmom_explicit_update,
1299                        &ymom_explicit_update,
1301                        &max_speed_array)) {
1302    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
1303    return NULL;
1304  }
1305  number_of_elements = stage_edge_values -> dimensions[0];
1306  call++;//a static local variable to which already_computed_flux is compared
1307  //set explicit_update to zero for all conserved_quantities.
1308  //This assumes compute_fluxes called before forcing terms
1309  for (k=0; k<number_of_elements; k++) {
1310    ((double *) stage_explicit_update -> data)[k]=0.0;
1311    ((double *) xmom_explicit_update -> data)[k]=0.0;
1312    ((double *) ymom_explicit_update -> data)[k]=0.0;
1313  }
1314
1315  //Loop through neighbours and compute edge flux for each
1316  for (k=0; k<number_of_elements; k++) {
1317    for (i=0; i<3; i++) {
1318      ki = k*3+i;
1319      if (((long *) already_computed_flux->data)[ki]==call)//we've already computed the flux across this edge
1320        continue;
1321      ql[0] = ((double *) stage_edge_values -> data)[ki];
1322      ql[1] = ((double *) xmom_edge_values -> data)[ki];
1323      ql[2] = ((double *) ymom_edge_values -> data)[ki];
1324      zl =    ((double *) bed_edge_values -> data)[ki];
1325
1326      //Quantities at neighbour on nearest face
1327      n = ((long *) neighbours -> data)[ki];
1328      if (n < 0) {
1329        m = -n-1; //Convert negative flag to index
1330        qr[0] = ((double *) stage_boundary_values -> data)[m];
1331        qr[1] = ((double *) xmom_boundary_values -> data)[m];
1332        qr[2] = ((double *) ymom_boundary_values -> data)[m];
1333        zr = zl; //Extend bed elevation to boundary
1334      } else {
1335        m = ((long *) neighbour_edges -> data)[ki];
1336        nm = n*3+m;
1337        qr[0] = ((double *) stage_edge_values -> data)[nm];
1338        qr[1] = ((double *) xmom_edge_values -> data)[nm];
1339        qr[2] = ((double *) ymom_edge_values -> data)[nm];
1340        zr =    ((double *) bed_edge_values -> data)[nm];
1341      }
1342      // Outward pointing normal vector
1343      // normal = domain.normals[k, 2*i:2*i+2]
1344      ki2 = 2*ki; //k*6 + i*2
1345      normal[0] = ((double *) normals -> data)[ki2];
1346      normal[1] = ((double *) normals -> data)[ki2+1];
1347      //Edge flux computation
1348      flux_function_central(ql, qr, zl, zr,
1349                    normal[0], normal[1],
1350                    epsilon, H0, g,
1351                    edgeflux, &max_speed);
1352      //update triangle k
1353      ((long *) already_computed_flux->data)[ki]=call;
1354      ((double *) stage_explicit_update -> data)[k] -= edgeflux[0]*((double *) edgelengths -> data)[ki];
1355      ((double *) xmom_explicit_update -> data)[k] -= edgeflux[1]*((double *) edgelengths -> data)[ki];
1356      ((double *) ymom_explicit_update -> data)[k] -= edgeflux[2]*((double *) edgelengths -> data)[ki];
1357      //update the neighbour n
1358      if (n>=0){
1359        ((long *) already_computed_flux->data)[nm]=call;
1360        ((double *) stage_explicit_update -> data)[n] += edgeflux[0]*((double *) edgelengths -> data)[nm];
1361        ((double *) xmom_explicit_update -> data)[n] += edgeflux[1]*((double *) edgelengths -> data)[nm];
1362        ((double *) ymom_explicit_update -> data)[n] += edgeflux[2]*((double *) edgelengths -> data)[nm];
1363      }
1364      ///for (j=0; j<3; j++) {
1365        ///flux[j] -= edgeflux[j]*((double *) edgelengths -> data)[ki];
1366        ///}
1367        //Update timestep
1368        //timestep = min(timestep, domain.radii[k]/max_speed)
1369        //FIXME: SR Add parameter for CFL condition
1370    if ( ((long *) tri_full_flag -> data)[k] == 1) {
1371            if (max_speed > epsilon) {
1372                timestep = min(timestep, ((double *) radii -> data)[k]/max_speed);
1373                //maxspeed in flux_function is calculated as max(|u+a|,|u-a|)
1374                if (n>=0)
1375                    timestep = min(timestep, ((double *) radii -> data)[n]/max_speed);
1376            }
1377    }
1378    } // end for i
1379
1380    //Normalise by area and store for when all conserved
1381    //quantities get updated
1382    ((double *) stage_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
1383    ((double *) xmom_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
1384    ((double *) ymom_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
1385
1386
1387    //Keep track of maximal speeds
1388    ((double *) max_speed_array -> data)[k] = max_speed;
1389
1390  } //end for k
1391  return Py_BuildValue("d", timestep);
1392}
1393
1394
1395PyObject *compute_fluxes_ext_kinetic(PyObject *self, PyObject *args) {
1396  /*Compute all fluxes and the timestep suitable for all volumes
1397    in domain.
1398
1399    Compute total flux for each conserved quantity using "flux_function_central"
1400
1401    Fluxes across each edge are scaled by edgelengths and summed up
1402    Resulting flux is then scaled by area and stored in
1403    explicit_update for each of the three conserved quantities
1404    stage, xmomentum and ymomentum
1405
1406    The maximal allowable speed computed by the flux_function for each volume
1407    is converted to a timestep that must not be exceeded. The minimum of
1408    those is computed as the next overall timestep.
1409
1410    Python call:
1411    domain.timestep = compute_fluxes(timestep,
1412                                     domain.epsilon,
1413                                     domain.H0,
1414                                     domain.g,
1415                                     domain.neighbours,
1416                                     domain.neighbour_edges,
1417                                     domain.normals,
1418                                     domain.edgelengths,
1420                                     domain.areas,
1421                                     Stage.edge_values,
1422                                     Xmom.edge_values,
1423                                     Ymom.edge_values,
1424                                     Bed.edge_values,
1425                                     Stage.boundary_values,
1426                                     Xmom.boundary_values,
1427                                     Ymom.boundary_values,
1428                                     Stage.explicit_update,
1429                                     Xmom.explicit_update,
1430                                     Ymom.explicit_update,
1432
1433
1434    Post conditions:
1435      domain.explicit_update is reset to computed flux values
1436      domain.timestep is set to the largest step satisfying all volumes.
1437
1438
1439  */
1440
1441
1442  PyArrayObject *neighbours, *neighbour_edges,
1443    *normals, *edgelengths, *radii, *areas,
1444    *tri_full_flag,
1445    *stage_edge_values,
1446    *xmom_edge_values,
1447    *ymom_edge_values,
1448    *bed_edge_values,
1449    *stage_boundary_values,
1450    *xmom_boundary_values,
1451    *ymom_boundary_values,
1452    *stage_explicit_update,
1453    *xmom_explicit_update,
1454    *ymom_explicit_update,
1455    *already_computed_flux;//tracks whether the flux across an edge has already been computed
1456
1457
1458  //Local variables
1459  double timestep, max_speed, epsilon, g, H0;
1460  double normal[2], ql[3], qr[3], zl, zr;
1461  double edgeflux[3]; //Work arrays for summing up fluxes
1462
1463  int number_of_elements, k, i, m, n;
1464  int ki, nm=0, ki2; //Index shorthands
1465  static long call=1;
1466
1467
1468  // Convert Python arguments to C
1469  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOO",
1470                        &timestep,
1471                        &epsilon,
1472                        &H0,
1473                        &g,
1474                        &neighbours,
1475                        &neighbour_edges,
1476                        &normals,
1477                        &edgelengths, &radii, &areas,
1478                        &tri_full_flag,
1479                        &stage_edge_values,
1480                        &xmom_edge_values,
1481                        &ymom_edge_values,
1482                        &bed_edge_values,
1483                        &stage_boundary_values,
1484                        &xmom_boundary_values,
1485                        &ymom_boundary_values,
1486                        &stage_explicit_update,
1487                        &xmom_explicit_update,
1488                        &ymom_explicit_update,
1490    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
1491    return NULL;
1492  }
1493  number_of_elements = stage_edge_values -> dimensions[0];
1494  call++;//a static local variable to which already_computed_flux is compared
1495  //set explicit_update to zero for all conserved_quantities.
1496  //This assumes compute_fluxes called before forcing terms
1497  for (k=0; k<number_of_elements; k++) {
1498    ((double *) stage_explicit_update -> data)[k]=0.0;
1499    ((double *) xmom_explicit_update -> data)[k]=0.0;
1500    ((double *) ymom_explicit_update -> data)[k]=0.0;
1501  }
1502  //Loop through neighbours and compute edge flux for each
1503  for (k=0; k<number_of_elements; k++) {
1504    for (i=0; i<3; i++) {
1505      ki = k*3+i;
1506      if (((long *) already_computed_flux->data)[ki]==call)//we've already computed the flux across this edge
1507        continue;
1508      ql[0] = ((double *) stage_edge_values -> data)[ki];
1509      ql[1] = ((double *) xmom_edge_values -> data)[ki];
1510      ql[2] = ((double *) ymom_edge_values -> data)[ki];
1511      zl =    ((double *) bed_edge_values -> data)[ki];
1512
1513      //Quantities at neighbour on nearest face
1514      n = ((long *) neighbours -> data)[ki];
1515      if (n < 0) {
1516        m = -n-1; //Convert negative flag to index
1517        qr[0] = ((double *) stage_boundary_values -> data)[m];
1518        qr[1] = ((double *) xmom_boundary_values -> data)[m];
1519        qr[2] = ((double *) ymom_boundary_values -> data)[m];
1520        zr = zl; //Extend bed elevation to boundary
1521      } else {
1522        m = ((long *) neighbour_edges -> data)[ki];
1523        nm = n*3+m;
1524        qr[0] = ((double *) stage_edge_values -> data)[nm];
1525        qr[1] = ((double *) xmom_edge_values -> data)[nm];
1526        qr[2] = ((double *) ymom_edge_values -> data)[nm];
1527        zr =    ((double *) bed_edge_values -> data)[nm];
1528      }
1529      // Outward pointing normal vector
1530      // normal = domain.normals[k, 2*i:2*i+2]
1531      ki2 = 2*ki; //k*6 + i*2
1532      normal[0] = ((double *) normals -> data)[ki2];
1533      normal[1] = ((double *) normals -> data)[ki2+1];
1534      //Edge flux computation
1535      flux_function_kinetic(ql, qr, zl, zr,
1536                    normal[0], normal[1],
1537                    epsilon, H0, g,
1538                    edgeflux, &max_speed);
1539      //update triangle k
1540      ((long *) already_computed_flux->data)[ki]=call;
1541      ((double *) stage_explicit_update -> data)[k] -= edgeflux[0]*((double *) edgelengths -> data)[ki];
1542      ((double *) xmom_explicit_update -> data)[k] -= edgeflux[1]*((double *) edgelengths -> data)[ki];
1543      ((double *) ymom_explicit_update -> data)[k] -= edgeflux[2]*((double *) edgelengths -> data)[ki];
1544      //update the neighbour n
1545      if (n>=0){
1546        ((long *) already_computed_flux->data)[nm]=call;
1547        ((double *) stage_explicit_update -> data)[n] += edgeflux[0]*((double *) edgelengths -> data)[nm];
1548        ((double *) xmom_explicit_update -> data)[n] += edgeflux[1]*((double *) edgelengths -> data)[nm];
1549        ((double *) ymom_explicit_update -> data)[n] += edgeflux[2]*((double *) edgelengths -> data)[nm];
1550      }
1551      ///for (j=0; j<3; j++) {
1552        ///flux[j] -= edgeflux[j]*((double *) edgelengths -> data)[ki];
1553        ///}
1554        //Update timestep
1555        //timestep = min(timestep, domain.radii[k]/max_speed)
1556        //FIXME: SR Add parameter for CFL condition
1557    if ( ((long *) tri_full_flag -> data)[k] == 1) {
1558            if (max_speed > epsilon) {
1559                timestep = min(timestep, ((double *) radii -> data)[k]/max_speed);
1560                //maxspeed in flux_function is calculated as max(|u+a|,|u-a|)
1561                if (n>=0)
1562                    timestep = min(timestep, ((double *) radii -> data)[n]/max_speed);
1563            }
1564    }
1565    } // end for i
1566    //Normalise by area and store for when all conserved
1567    //quantities get updated
1568    ((double *) stage_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
1569    ((double *) xmom_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
1570    ((double *) ymom_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
1571  } //end for k
1572  return Py_BuildValue("d", timestep);
1573}
1574
1575PyObject *protect(PyObject *self, PyObject *args) {
1576  //
1577  //    protect(minimum_allowed_height, maximum_allowed_speed, wc, zc, xmomc, ymomc)
1578
1579
1580  PyArrayObject
1581  *wc,            //Stage at centroids
1582  *zc,            //Elevation at centroids
1583  *xmomc,         //Momentums at centroids
1584  *ymomc;
1585
1586
1587  int N;
1588  double minimum_allowed_height, maximum_allowed_speed, epsilon;
1589
1590  // Convert Python arguments to C
1591  if (!PyArg_ParseTuple(args, "dddOOOO",
1592                        &minimum_allowed_height,
1593                        &maximum_allowed_speed,
1594                        &epsilon,
1595                        &wc, &zc, &xmomc, &ymomc)) {
1596    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: protect could not parse input arguments");
1597    return NULL;
1598  }
1599
1600  N = wc -> dimensions[0];
1601
1602  _protect(N,
1603           minimum_allowed_height,
1604           maximum_allowed_speed,
1605           epsilon,
1606           (double*) wc -> data,
1607           (double*) zc -> data,
1608           (double*) xmomc -> data,
1609           (double*) ymomc -> data);
1610
1611  return Py_BuildValue("");
1612}
1613
1614
1615
1616PyObject *balance_deep_and_shallow(PyObject *self, PyObject *args) {
1617  //
1618  //    balance_deep_and_shallow(wc, zc, hc, wv, zv, hv,
1619  //                             xmomc, ymomc, xmomv, ymomv)
1620
1621
1622  PyArrayObject
1623    *wc,            //Stage at centroids
1624    *zc,            //Elevation at centroids
1625    *hc,            //Height at centroids
1626    *wv,            //Stage at vertices
1627    *zv,            //Elevation at vertices
1628    *hv,            //Depths at vertices
1629    *hvbar,         //h-Limited depths at vertices
1630    *xmomc,         //Momentums at centroids and vertices
1631    *ymomc,
1632    *xmomv,
1633    *ymomv;
1634
1635  PyObject *domain, *Tmp;
1636
1637  double alpha_balance = 2.0;
1638  double H0;
1639
1640  int N, limit2007; //, err;
1641
1642  // Convert Python arguments to C
1643  if (!PyArg_ParseTuple(args, "OOOOOOOOOOOO",
1644                        &domain,
1645                        &wc, &zc, &hc,
1646                        &wv, &zv, &hv, &hvbar,
1647                        &xmomc, &ymomc, &xmomv, &ymomv)) {
1648    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not parse input arguments");
1649    return NULL;
1650  }
1651
1652  // Pull out parameters
1653  Tmp = PyObject_GetAttrString(domain, "alpha_balance");
1654  if (!Tmp) {
1655    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not obtain object alpha_balance from domain");
1656    return NULL;
1657  }
1658  alpha_balance = PyFloat_AsDouble(Tmp);
1659  Py_DECREF(Tmp);
1660
1661
1662  Tmp = PyObject_GetAttrString(domain, "H0");
1663  if (!Tmp) {
1664    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not obtain object H0 from domain");
1665    return NULL;
1666  }
1667  H0 = PyFloat_AsDouble(Tmp);
1668  Py_DECREF(Tmp);
1669
1670
1671  Tmp = PyObject_GetAttrString(domain, "limit2007");
1672  if (!Tmp) {
1673    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not obtain object limit2007 from domain");
1674    return NULL;
1675  }
1676  limit2007 = PyInt_AsLong(Tmp);
1677  Py_DECREF(Tmp);
1678
1679
1680
1681
1682  N = wc -> dimensions[0];
1683
1684  _balance_deep_and_shallow(N,
1685                            (double*) wc -> data,
1686                            (double*) zc -> data,
1687                            (double*) hc -> data,
1688                            (double*) wv -> data,
1689                            (double*) zv -> data,
1690                            (double*) hv -> data,
1691                            (double*) hvbar -> data,
1692                            (double*) xmomc -> data,
1693                            (double*) ymomc -> data,
1694                            (double*) xmomv -> data,
1695                            (double*) ymomv -> data,
1696                            H0,
1697                            (int) limit2007,
1698                            alpha_balance);
1699
1700
1701  return Py_BuildValue("");
1702}
1703
1704
1705
1706PyObject *h_limiter(PyObject *self, PyObject *args) {
1707
1708  PyObject *domain, *Tmp;
1709  PyArrayObject
1710    *hv, *hc, //Depth at vertices and centroids
1711    *hvbar,   //Limited depth at vertices (return values)
1712    *neighbours;
1713
1714  int k, i, n, N, k3;
1715  int dimensions[2];
1716  double beta_h; //Safety factor (see config.py)
1717  double *hmin, *hmax, hn;
1718
1719  // Convert Python arguments to C
1720  if (!PyArg_ParseTuple(args, "OOO", &domain, &hc, &hv)) {
1721    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: h_limiter could not parse input arguments");
1722    return NULL;
1723  }
1724
1725  neighbours = get_consecutive_array(domain, "neighbours");
1726
1727  //Get safety factor beta_h
1728  Tmp = PyObject_GetAttrString(domain, "beta_h");
1729  if (!Tmp) {
1730    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: h_limiter could not obtain object beta_h from domain");
1731    return NULL;
1732  }
1733  beta_h = PyFloat_AsDouble(Tmp);
1734
1735  Py_DECREF(Tmp);
1736
1737  N = hc -> dimensions[0];
1738
1739  //Create hvbar
1740  dimensions[0] = N;
1741  dimensions[1] = 3;
1742  hvbar = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
1743
1744
1745  //Find min and max of this and neighbour's centroid values
1746  hmin = malloc(N * sizeof(double));
1747  hmax = malloc(N * sizeof(double));
1748  for (k=0; k<N; k++) {
1749    k3=k*3;
1750
1751    hmin[k] = ((double*) hc -> data)[k];
1752    hmax[k] = hmin[k];
1753
1754    for (i=0; i<3; i++) {
1755      n = ((long*) neighbours -> data)[k3+i];
1756
1757      //Initialise hvbar with values from hv
1758      ((double*) hvbar -> data)[k3+i] = ((double*) hv -> data)[k3+i];
1759
1760      if (n >= 0) {
1761        hn = ((double*) hc -> data)[n]; //Neighbour's centroid value
1762
1763        hmin[k] = min(hmin[k], hn);
1764        hmax[k] = max(hmax[k], hn);
1765      }
1766    }
1767  }
1768
1769  // Call underlying standard routine
1770  _limit(N, beta_h, (double*) hc -> data, (double*) hvbar -> data, hmin, hmax);
1771
1772  // // //Py_DECREF(domain); //FIXME: NEcessary?
1773  free(hmin);
1774  free(hmax);
1775
1776  //return result using PyArray to avoid memory leak
1777  return PyArray_Return(hvbar);
1778  //return Py_BuildValue("");
1779}
1780
1781PyObject *h_limiter_sw(PyObject *self, PyObject *args) {
1782  //a faster version of h_limiter above
1783  PyObject *domain, *Tmp;
1784  PyArrayObject
1785    *hv, *hc, //Depth at vertices and centroids
1786    *hvbar,   //Limited depth at vertices (return values)
1787    *neighbours;
1788
1789  int k, i, N, k3,k0,k1,k2;
1790  int dimensions[2];
1791  double beta_h; //Safety factor (see config.py)
1792  double hmin, hmax, dh[3];
1793// Convert Python arguments to C
1794  if (!PyArg_ParseTuple(args, "OOO", &domain, &hc, &hv)) {
1795    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: h_limiter_sw could not parse input arguments");
1796    return NULL;
1797  }
1798  neighbours = get_consecutive_array(domain, "neighbours");
1799
1800  //Get safety factor beta_h
1801  Tmp = PyObject_GetAttrString(domain, "beta_h");
1802  if (!Tmp) {
1803    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: h_limiter_sw could not obtain object beta_h from domain");
1804    return NULL;
1805  }
1806  beta_h = PyFloat_AsDouble(Tmp);
1807
1808  Py_DECREF(Tmp);
1809
1810  N = hc -> dimensions[0];
1811
1812  //Create hvbar
1813  dimensions[0] = N;
1814  dimensions[1] = 3;
1815  hvbar = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
1816  for (k=0;k<N;k++){
1817    k3=k*3;
1818    //get the ids of the neighbours
1819    k0 = ((long*) neighbours -> data)[k3];
1820    k1 = ((long*) neighbours -> data)[k3+1];
1821    k2 = ((long*) neighbours -> data)[k3+2];
1822    //set hvbar provisionally
1823    for (i=0;i<3;i++){
1824      ((double*) hvbar -> data)[k3+i] = ((double*) hv -> data)[k3+i];
1825      dh[i]=((double*) hvbar -> data)[k3+i]-((double*) hc -> data)[k];
1826    }
1827    hmin=((double*) hc -> data)[k];
1828    hmax=hmin;
1829    if (k0>=0){
1830      hmin=min(hmin,((double*) hc -> data)[k0]);
1831      hmax=max(hmax,((double*) hc -> data)[k0]);
1832    }
1833    if (k1>=0){
1834      hmin=min(hmin,((double*) hc -> data)[k1]);
1835      hmax=max(hmax,((double*) hc -> data)[k1]);
1836    }
1837    if (k2>=0){
1838      hmin=min(hmin,((double*) hc -> data)[k2]);
1839      hmax=max(hmax,((double*) hc -> data)[k2]);
1840    }
1841    hmin-=((double*) hc -> data)[k];
1842    hmax-=((double*) hc -> data)[k];
1844    for (i=0;i<3;i++)
1845      ((double*) hvbar -> data)[k3+i] = ((double*) hc -> data)[k]+dh[i];
1846  }
1847  return PyArray_Return(hvbar);
1848}
1849
1850PyObject *assign_windfield_values(PyObject *self, PyObject *args) {
1851  //
1852  //      assign_windfield_values(xmom_update, ymom_update,
1853  //                              s_vec, phi_vec, self.const)
1854
1855
1856
1857  PyArrayObject   //(one element per triangle)
1858  *s_vec,         //Speeds
1859  *phi_vec,       //Bearings
1860  *xmom_update,   //Momentum updates
1861  *ymom_update;
1862
1863
1864  int N;
1865  double cw;
1866
1867  // Convert Python arguments to C
1868  if (!PyArg_ParseTuple(args, "OOOOd",
1869                        &xmom_update,
1870                        &ymom_update,
1871                        &s_vec, &phi_vec,
1872                        &cw)) {
1873    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: assign_windfield_values could not parse input arguments");
1874    return NULL;
1875  }
1876
1877
1878  N = xmom_update -> dimensions[0];
1879
1880  _assign_wind_field_values(N,
1881           (double*) xmom_update -> data,
1882           (double*) ymom_update -> data,
1883           (double*) s_vec -> data,
1884           (double*) phi_vec -> data,
1885           cw);
1886
1887  return Py_BuildValue("");
1888}
1889
1890
1891
1892
1893//////////////////////////////////////////
1894// Method table for python module
1895static struct PyMethodDef MethodTable[] = {
1896  /* The cast of the function is necessary since PyCFunction values
1897   * only take two PyObject* parameters, and rotate() takes
1898   * three.
1899   */
1900
1901  {"rotate", (PyCFunction)rotate, METH_VARARGS | METH_KEYWORDS, "Print out"},
1902  {"extrapolate_second_order_sw", extrapolate_second_order_sw, METH_VARARGS, "Print out"},
1903  {"compute_fluxes_ext_central", compute_fluxes_ext_central, METH_VARARGS, "Print out"},
1904  {"compute_fluxes_ext_kinetic", compute_fluxes_ext_kinetic, METH_VARARGS, "Print out"},
1905  {"gravity", gravity, METH_VARARGS, "Print out"},
1906  {"manning_friction", manning_friction, METH_VARARGS, "Print out"},
1907  {"balance_deep_and_shallow", balance_deep_and_shallow,
1908   METH_VARARGS, "Print out"},
1909  {"h_limiter", h_limiter,
1910   METH_VARARGS, "Print out"},
1911  {"h_limiter_sw", h_limiter_sw,
1912   METH_VARARGS, "Print out"},
1913  {"protect", protect, METH_VARARGS | METH_KEYWORDS, "Print out"},
1914  {"assign_windfield_values", assign_windfield_values,
1915   METH_VARARGS | METH_KEYWORDS, "Print out"},
1916  //{"distribute_to_vertices_and_edges",
1917  // distribute_to_vertices_and_edges, METH_VARARGS},
1918  //{"update_conserved_quantities",
1919  // update_conserved_quantities, METH_VARARGS},
1920  //{"set_initialcondition",
1921  // set_initialcondition, METH_VARARGS},
1922  {NULL, NULL}
1923};
1924
1925// Module initialisation
1926void initshallow_water_ext(void){
1927  Py_InitModule("shallow_water_ext", MethodTable);
1928
1929  import_array();     //Necessary for handling of NumPY structures
1930}
Note: See TracBrowser for help on using the repository browser.