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

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

Added set_maxmum_allowed_speed to zero to simulate behaviour of ANUGA as it was in August 2005. This provides the exact same look and feel as was the case when we first did the Okushiri validation.

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