source: anuga_core/source/anuga/pyvolution/shallow_water_ext.c @ 3546

Last change on this file since 3546 was 3246, checked in by jack, 19 years ago

Added skeleton SCons build files.

File size: 44.9 KB
Line 
1// Python - C extension module for shallow_water.py
2//
3// To compile (Python2.3):
4//  gcc -c domain_ext.c -I/usr/include/python2.3 -o domain_ext.o -Wall -O
5//  gcc -shared domain_ext.o  -o domain_ext.so
6//
7// or use python compile.py
8//
9// See the module shallow_water.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(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
235void _manning_friction(double g, double eps, int N,
236                       double* w, double* z,
237                       double* uh, double* vh,
238                       double* eta, double* xmom, double* ymom) {
239
240  int k;
241  double S, h;
242
243  for (k=0; k<N; k++) {
244    if (eta[k] > eps) {
245      h = w[k]-z[k];
246      if (h >= eps) {
247        S = -g * eta[k]*eta[k] * sqrt((uh[k]*uh[k] + vh[k]*vh[k]));
248        S /= pow(h, 7.0/3);      //Expensive (on Ole's home computer)
249        //S /= exp(7.0/3.0*log(h));      //seems to save about 15% over manning_friction
250        //S /= h*h*(1 + h/3.0 - h*h/9.0); //FIXME: Could use a Taylor expansion
251
252
253        //Update momentum
254        xmom[k] += S*uh[k];
255        ymom[k] += S*vh[k];
256      }
257    }
258  }
259}
260
261
262/*
263void _manning_friction_explicit(double g, double eps, int N,
264                       double* w, double* z,
265                       double* uh, double* vh,
266                       double* eta, double* xmom, double* ymom) {
267
268  int k;
269  double S, h;
270
271  for (k=0; k<N; k++) {
272    if (eta[k] > eps) {
273      h = w[k]-z[k];
274      if (h >= eps) {
275        S = -g * eta[k]*eta[k] * sqrt((uh[k]*uh[k] + vh[k]*vh[k]));
276        S /= pow(h, 7.0/3);      //Expensive (on Ole's home computer)
277        //S /= exp(7.0/3.0*log(h));      //seems to save about 15% over manning_friction
278        //S /= h*h*(1 + h/3.0 - h*h/9.0); //FIXME: Could use a Taylor expansion
279
280
281        //Update momentum
282        xmom[k] += S*uh[k];
283        ymom[k] += S*vh[k];
284      }
285    }
286  }
287}
288*/
289
290int _balance_deep_and_shallow(int N,
291                              double* wc,
292                              double* zc,
293                              double* hc,
294                              double* wv,
295                              double* zv,
296                              double* hv,
297                              double* hvbar,
298                              double* xmomc,
299                              double* ymomc,
300                              double* xmomv,
301                              double* ymomv) {
302
303  int k, k3, i;
304  double dz, hmin, alpha;
305
306  //Compute linear combination between w-limited stages and
307  //h-limited stages close to the bed elevation.
308
309  for (k=0; k<N; k++) {
310    // Compute maximal variation in bed elevation
311    // This quantitiy is
312    //     dz = max_i abs(z_i - z_c)
313    // and it is independent of dimension
314    // In the 1d case zc = (z0+z1)/2
315    // In the 2d case zc = (z0+z1+z2)/3
316
317    k3 = 3*k;
318
319    //FIXME: Try with this one precomputed
320    dz = 0.0;
321    hmin = hv[k3];
322    for (i=0; i<3; i++) {
323      dz = max(dz, fabs(zv[k3+i]-zc[k]));
324      hmin = min(hmin, hv[k3+i]);
325    }
326
327
328    //Create alpha in [0,1], where alpha==0 means using the h-limited
329    //stage and alpha==1 means using the w-limited stage as
330    //computed by the gradient limiter (both 1st or 2nd order)
331    //
332    //If hmin > dz/2 then alpha = 1 and the bed will have no effect
333    //If hmin < 0 then alpha = 0 reverting to constant height above bed.
334
335
336    if (dz > 0.0)
337      //if (hmin<0.0)
338      //        alpha = 0.0;
339      //else
340      //  alpha = max( min( hc[k]/dz, 1.0), 0.0 );
341      alpha = max( min( 2.0*hmin/dz, 1.0), 0.0 );
342    else
343      alpha = 1.0;  //Flat bed
344
345
346    //printf("dz = %.3f, alpha = %.8f\n", dz, alpha);
347
348    //  Let
349    //
350    //    wvi be the w-limited stage (wvi = zvi + hvi)
351    //    wvi- be the h-limited state (wvi- = zvi + hvi-)
352    //
353    //
354    //  where i=0,1,2 denotes the vertex ids
355    //
356    //  Weighted balance between w-limited and h-limited stage is
357    //
358    //    wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi)
359    //
360    //  It follows that the updated wvi is
361    //    wvi := zvi + (1-alpha)*hvi- + alpha*hvi
362    //
363    //   Momentum is balanced between constant and limited
364
365    if (alpha < 1) {
366      for (i=0; i<3; i++) {
367         wv[k3+i] = zv[k3+i] + (1-alpha)*hvbar[k3+i] + alpha*hv[k3+i];
368
369        //Update momentum as a linear combination of
370        //xmomc and ymomc (shallow) and momentum
371        //from extrapolator xmomv and ymomv (deep).
372        xmomv[k3+i] = (1-alpha)*xmomc[k] + alpha*xmomv[k3+i];
373        ymomv[k3+i] = (1-alpha)*ymomc[k] + alpha*ymomv[k3+i];
374      }
375    }
376  }
377  return 0;
378}
379
380
381
382int _protect(int N,
383             double minimum_allowed_height,
384             double maximum_allowed_speed,
385             double epsilon,
386             double* wc,
387             double* zc,
388             double* xmomc,
389             double* ymomc) {
390
391  int k;
392  double hc;
393  double u, v, reduced_speed;
394
395  //Protect against initesimal and negative heights
396  for (k=0; k<N; k++) {
397    hc = wc[k] - zc[k];
398
399    if (hc < minimum_allowed_height) {
400
401      //Old code: Set momentum to zero and ensure h is non negative
402      //xmomc[k] = 0.0;
403      //ymomc[k] = 0.0;
404      //if (hc <= 0.0) wc[k] = zc[k];
405
406
407      //New code: Adjust momentum to guarantee speeds are physical
408      //          ensure h is non negative
409      //FIXME (Ole): This is only implemented in this C extension and
410      //             has no Python equivalent
411      if (hc <= 0.0) {
412        wc[k] = zc[k];
413        xmomc[k] = 0.0;
414        ymomc[k] = 0.0;
415      } else {
416        //Reduce excessive speeds derived from division by small hc
417        u = xmomc[k]/hc;
418        if (fabs(u) > maximum_allowed_speed) {
419          reduced_speed = maximum_allowed_speed * u/fabs(u);
420          //printf("Speed (u) has been reduced from %.3f to %.3f\n",
421          //     u, reduced_speed);
422          xmomc[k] = reduced_speed * hc;
423        }
424
425        v = ymomc[k]/hc;
426        if (fabs(v) > maximum_allowed_speed) {
427          reduced_speed = maximum_allowed_speed * v/fabs(v);
428          //printf("Speed (v) has been reduced from %.3f to %.3f\n",
429          //     v, reduced_speed);
430          ymomc[k] = reduced_speed * hc;
431        }
432      }
433    }
434  }
435  return 0;
436}
437
438
439
440
441int _assign_wind_field_values(int N,
442                              double* xmom_update,
443                              double* ymom_update,
444                              double* s_vec,
445                              double* phi_vec,
446                              double cw) {
447
448  //Assign windfield values to momentum updates
449
450  int k;
451  double S, s, phi, u, v;
452
453  for (k=0; k<N; k++) {
454
455    s = s_vec[k];
456    phi = phi_vec[k];
457
458    //Convert to radians
459    phi = phi*pi/180;
460
461    //Compute velocity vector (u, v)
462    u = s*cos(phi);
463    v = s*sin(phi);
464
465    //Compute wind stress
466    S = cw * sqrt(u*u + v*v);
467    xmom_update[k] += S*u;
468    ymom_update[k] += S*v;
469  }
470  return 0;
471}
472
473
474
475///////////////////////////////////////////////////////////////////
476// Gateways to Python
477
478PyObject *gravity(PyObject *self, PyObject *args) {
479  //
480  //  gravity(g, h, v, x, xmom, ymom)
481  //
482
483
484  PyArrayObject *h, *v, *x, *xmom, *ymom;
485  int k, i, N, k3, k6;
486  double g, avg_h, zx, zy;
487  double x0, y0, x1, y1, x2, y2, z0, z1, z2;
488
489  if (!PyArg_ParseTuple(args, "dOOOOO",
490                        &g, &h, &v, &x,
491                        &xmom, &ymom))
492    return NULL;
493
494  N = h -> dimensions[0];
495  for (k=0; k<N; k++) {
496    k3 = 3*k;  // base index
497    k6 = 6*k;  // base index
498
499    avg_h = 0.0;
500    for (i=0; i<3; i++) {
501      avg_h += ((double *) h -> data)[k3+i];
502    }
503    avg_h /= 3;
504
505
506    //Compute bed slope
507    x0 = ((double*) x -> data)[k6 + 0];
508    y0 = ((double*) x -> data)[k6 + 1];
509    x1 = ((double*) x -> data)[k6 + 2];
510    y1 = ((double*) x -> data)[k6 + 3];
511    x2 = ((double*) x -> data)[k6 + 4];
512    y2 = ((double*) x -> data)[k6 + 5];
513
514
515    z0 = ((double*) v -> data)[k3 + 0];
516    z1 = ((double*) v -> data)[k3 + 1];
517    z2 = ((double*) v -> data)[k3 + 2];
518
519    _gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2, &zx, &zy);
520
521    //Update momentum
522    ((double*) xmom -> data)[k] += -g*zx*avg_h;
523    ((double*) ymom -> data)[k] += -g*zy*avg_h;
524  }
525
526  return Py_BuildValue("");
527}
528
529
530PyObject *manning_friction(PyObject *self, PyObject *args) {
531  //
532  // manning_friction(g, eps, h, uh, vh, eta, xmom_update, ymom_update)
533  //
534
535
536  PyArrayObject *w, *z, *uh, *vh, *eta, *xmom, *ymom;
537  int N;
538  double g, eps;
539
540  if (!PyArg_ParseTuple(args, "ddOOOOOOO",
541                        &g, &eps, &w, &z, &uh, &vh, &eta,
542                        &xmom, &ymom))
543    return NULL;
544
545  N = w -> dimensions[0];
546  _manning_friction(g, eps, N,
547                    (double*) w -> data,
548                    (double*) z -> data,
549                    (double*) uh -> data,
550                    (double*) vh -> data,
551                    (double*) eta -> data,
552                    (double*) xmom -> data,
553                    (double*) ymom -> data);
554
555  return Py_BuildValue("");
556}
557
558
559/*
560PyObject *manning_friction_explicit(PyObject *self, PyObject *args) {
561  //
562  // manning_friction_explicit(g, eps, h, uh, vh, eta, xmom_update, ymom_update)
563  //
564
565
566  PyArrayObject *w, *z, *uh, *vh, *eta, *xmom, *ymom;
567  int N;
568  double g, eps;
569
570  if (!PyArg_ParseTuple(args, "ddOOOOOOO",
571                        &g, &eps, &w, &z, &uh, &vh, &eta,
572                        &xmom, &ymom))
573    return NULL;
574
575  N = w -> dimensions[0];
576  _manning_friction_explicit(g, eps, N,
577                    (double*) w -> data,
578                    (double*) z -> data,
579                    (double*) uh -> data,
580                    (double*) vh -> data,
581                    (double*) eta -> data,
582                    (double*) xmom -> data,
583                    (double*) ymom -> data);
584
585  return Py_BuildValue("");
586}
587*/
588
589PyObject *extrapolate_second_order_sw(PyObject *self, PyObject *args) {
590  /*Compute the vertex values based on a linear reconstruction on each triangle
591    These values are calculated as follows:
592    1) For each triangle not adjacent to a boundary, we consider the auxiliary triangle
593    formed by the centroids of its three neighbours.
594    2) For each conserved quantity, we integrate around the auxiliary triangle's boundary the product
595    of the quantity and the outward normal vector. Dividing by the triangle area gives (a,b), the average
596    of the vector (q_x,q_y) on the auxiliary triangle. We suppose that the linear reconstruction on the
597    original triangle has gradient (a,b).
598    3) Provisional vertex junmps dqv[0,1,2] are computed and these are then limited by calling the functions
599    find_qmin_and_qmax and limit_gradient
600
601    Python call:
602    extrapolate_second_order_sw(domain.surrogate_neighbours,
603                                domain.number_of_boundaries
604                                domain.centroid_coordinates,
605                                Stage.centroid_values
606                                Xmom.centroid_values
607                                Ymom.centroid_values
608                                domain.vertex_coordinates,
609                                Stage.vertex_values,
610                                Xmom.vertex_values,
611                                Ymom.vertex_values)
612
613    Post conditions:
614            The vertices of each triangle have values from a limited linear reconstruction
615            based on centroid values
616
617  */
618  PyArrayObject *surrogate_neighbours,
619    *number_of_boundaries,
620    *centroid_coordinates,
621    *stage_centroid_values,
622    *xmom_centroid_values,
623    *ymom_centroid_values,
624    *vertex_coordinates,
625    *stage_vertex_values,
626    *xmom_vertex_values,
627    *ymom_vertex_values;
628  PyObject *domain, *Tmp;
629  //Local variables
630  double a, b;//gradient vector, not stored but used to calculate vertex values from centroids
631  int number_of_elements,k,k0,k1,k2,k3,k6,coord_index,i;
632  double x,y,x0,y0,x1,y1,x2,y2,xv0,yv0,xv1,yv1,xv2,yv2;//vertices of the auxiliary triangle
633  double dx1,dx2,dy1,dy2,dxv0,dxv1,dxv2,dyv0,dyv1,dyv2,dq0,dq1,dq2,area2;
634  double dqv[3], qmin, qmax, beta_w;//provisional jumps from centroids to v'tices and safety factor re limiting
635  //by which these jumps are limited
636  // Convert Python arguments to C
637  if (!PyArg_ParseTuple(args, "OOOOOOOOOOO",
638                        &domain,
639                        &surrogate_neighbours,
640                        &number_of_boundaries,
641                        &centroid_coordinates,
642                        &stage_centroid_values,
643                        &xmom_centroid_values,
644                        &ymom_centroid_values,
645                        &vertex_coordinates,
646                        &stage_vertex_values,
647                        &xmom_vertex_values,
648                        &ymom_vertex_values)) {
649    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
650    return NULL;
651  }
652
653  //get the safety factor beta_w, set in the config.py file. This is used in the limiting process
654  Tmp = PyObject_GetAttrString(domain, "beta_w");
655  if (!Tmp)
656    return NULL;
657  beta_w = PyFloat_AsDouble(Tmp);
658  Py_DECREF(Tmp);
659  number_of_elements = stage_centroid_values -> dimensions[0];
660  for (k=0; k<number_of_elements; k++) {
661    k3=k*3;
662    k6=k*6;
663
664    if (((long *) number_of_boundaries->data)[k]==3){/*no neighbours, set gradient on the triangle to zero*/
665      ((double *) stage_vertex_values->data)[k3]=((double *)stage_centroid_values->data)[k];
666      ((double *) stage_vertex_values->data)[k3+1]=((double *)stage_centroid_values->data)[k];
667      ((double *) stage_vertex_values->data)[k3+2]=((double *)stage_centroid_values->data)[k];
668      ((double *) xmom_vertex_values->data)[k3]=((double *)xmom_centroid_values->data)[k];
669      ((double *) xmom_vertex_values->data)[k3+1]=((double *)xmom_centroid_values->data)[k];
670      ((double *) xmom_vertex_values->data)[k3+2]=((double *)xmom_centroid_values->data)[k];
671      ((double *) ymom_vertex_values->data)[k3]=((double *)ymom_centroid_values->data)[k];
672      ((double *) ymom_vertex_values->data)[k3+1]=((double *)ymom_centroid_values->data)[k];
673      ((double *) ymom_vertex_values->data)[k3+2]=((double *)ymom_centroid_values->data)[k];
674      continue;
675    }
676    else {//we will need centroid coordinates and vertex coordinates of the triangle
677      //get the vertex coordinates of the FV triangle
678      xv0=((double *)vertex_coordinates->data)[k6]; yv0=((double *)vertex_coordinates->data)[k6+1];
679      xv1=((double *)vertex_coordinates->data)[k6+2]; yv1=((double *)vertex_coordinates->data)[k6+3];
680      xv2=((double *)vertex_coordinates->data)[k6+4]; yv2=((double *)vertex_coordinates->data)[k6+5];
681      //get the centroid coordinates of the FV triangle
682      coord_index=2*k;
683      x=((double *)centroid_coordinates->data)[coord_index];
684      y=((double *)centroid_coordinates->data)[coord_index+1];
685      //store x- and y- differentials for the vertices of the FV triangle relative to the centroid
686      dxv0=xv0-x; dxv1=xv1-x; dxv2=xv2-x;
687      dyv0=yv0-y; dyv1=yv1-y; dyv2=yv2-y;
688    }
689    if (((long *)number_of_boundaries->data)[k]<=1){
690      //if no boundaries, auxiliary triangle is formed from the centroids of the three neighbours
691      //if one boundary, auxiliary triangle is formed from this centroid and its two neighbours
692      k0=((long *)surrogate_neighbours->data)[k3];
693      k1=((long *)surrogate_neighbours->data)[k3+1];
694      k2=((long *)surrogate_neighbours->data)[k3+2];
695      //get the auxiliary triangle's vertex coordinates (really the centroids of neighbouring triangles)
696      coord_index=2*k0;
697      x0=((double *)centroid_coordinates->data)[coord_index];
698      y0=((double *)centroid_coordinates->data)[coord_index+1];
699      coord_index=2*k1;
700      x1=((double *)centroid_coordinates->data)[coord_index];
701      y1=((double *)centroid_coordinates->data)[coord_index+1];
702      coord_index=2*k2;
703      x2=((double *)centroid_coordinates->data)[coord_index];
704      y2=((double *)centroid_coordinates->data)[coord_index+1];
705      //store x- and y- differentials for the vertices of the auxiliary triangle
706      dx1=x1-x0; dx2=x2-x0;
707      dy1=y1-y0; dy2=y2-y0;
708      //calculate 2*area of the auxiliary triangle
709      area2 = dy2*dx1 - dy1*dx2;//the triangle is guaranteed to be counter-clockwise
710      //If the mesh is 'weird' near the boundary, the trianlge might be flat or clockwise:
711      if (area2<=0)
712        return NULL;
713
714      //### stage ###
715      //calculate the difference between vertex 0 of the auxiliary triangle and the FV triangle centroid
716      dq0=((double *)stage_centroid_values->data)[k0]-((double *)stage_centroid_values->data)[k];
717      //calculate differentials between the vertices of the auxiliary triangle
718      dq1=((double *)stage_centroid_values->data)[k1]-((double *)stage_centroid_values->data)[k0];
719      dq2=((double *)stage_centroid_values->data)[k2]-((double *)stage_centroid_values->data)[k0];
720      //calculate the gradient of stage on the auxiliary triangle
721      a = dy2*dq1 - dy1*dq2;
722      a /= area2;
723      b = dx1*dq2 - dx2*dq1;
724      b /= area2;
725      //calculate provisional jumps in stage from the centroid of the FV tri to its vertices, to be limited
726      dqv[0]=a*dxv0+b*dyv0;
727      dqv[1]=a*dxv1+b*dyv1;
728      dqv[2]=a*dxv2+b*dyv2;
729      //now we want to find min and max of the centroid and the vertices of the auxiliary triangle
730      //and compute jumps from the centroid to the min and max
731      find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
732      limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
733      for (i=0;i<3;i++)
734        ((double *)stage_vertex_values->data)[k3+i]=((double *)stage_centroid_values->data)[k]+dqv[i];
735
736      //### xmom ###
737      //calculate the difference between vertex 0 of the auxiliary triangle and the FV triangle centroid
738      dq0=((double *)xmom_centroid_values->data)[k0]-((double *)xmom_centroid_values->data)[k];
739      //calculate differentials between the vertices of the auxiliary triangle
740      dq1=((double *)xmom_centroid_values->data)[k1]-((double *)xmom_centroid_values->data)[k0];
741      dq2=((double *)xmom_centroid_values->data)[k2]-((double *)xmom_centroid_values->data)[k0];
742      //calculate the gradient of xmom on the auxiliary triangle
743      a = dy2*dq1 - dy1*dq2;
744      a /= area2;
745      b = dx1*dq2 - dx2*dq1;
746      b /= area2;
747      //calculate provisional jumps in stage from the centroid of the FV tri to its vertices, to be limited
748      dqv[0]=a*dxv0+b*dyv0;
749      dqv[1]=a*dxv1+b*dyv1;
750      dqv[2]=a*dxv2+b*dyv2;
751      //now we want to find min and max of the centroid and the vertices of the auxiliary triangle
752      //and compute jumps from the centroid to the min and max
753      find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
754      limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
755      for (i=0;i<3;i++)
756        ((double *)xmom_vertex_values->data)[k3+i]=((double *)xmom_centroid_values->data)[k]+dqv[i];
757
758      //### ymom ###
759      //calculate the difference between vertex 0 of the auxiliary triangle and the FV triangle centroid
760      dq0=((double *)ymom_centroid_values->data)[k0]-((double *)ymom_centroid_values->data)[k];
761      //calculate differentials between the vertices of the auxiliary triangle
762      dq1=((double *)ymom_centroid_values->data)[k1]-((double *)ymom_centroid_values->data)[k0];
763      dq2=((double *)ymom_centroid_values->data)[k2]-((double *)ymom_centroid_values->data)[k0];
764      //calculate the gradient of xmom on the auxiliary triangle
765      a = dy2*dq1 - dy1*dq2;
766      a /= area2;
767      b = dx1*dq2 - dx2*dq1;
768      b /= area2;
769      //calculate provisional jumps in stage from the centroid of the FV tri to its vertices, to be limited
770      dqv[0]=a*dxv0+b*dyv0;
771      dqv[1]=a*dxv1+b*dyv1;
772      dqv[2]=a*dxv2+b*dyv2;
773      //now we want to find min and max of the centroid and the vertices of the auxiliary triangle
774      //and compute jumps from the centroid to the min and max
775      find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
776      limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
777      for (i=0;i<3;i++)
778        ((double *)ymom_vertex_values->data)[k3+i]=((double *)ymom_centroid_values->data)[k]+dqv[i];
779    }//if (number_of_boundaries[k]<=1)
780    else{//number_of_boundaries==2
781      //one internal neighbour and gradient is in direction of the neighbour's centroid
782      //find the only internal neighbour
783      for (k2=k3;k2<k3+3;k2++){//k2 just indexes the edges of triangle k
784        if (((long *)surrogate_neighbours->data)[k2]!=k)//find internal neighbour of triabngle k
785          break;
786      }
787      if ((k2==k3+3))//if we didn't find an internal neighbour
788        return NULL;//error
789      k1=((long *)surrogate_neighbours->data)[k2];
790      //the coordinates of the triangle are already (x,y). Get centroid of the neighbour (x1,y1)
791      coord_index=2*k1;
792      x1=((double *)centroid_coordinates->data)[coord_index];
793      y1=((double *)centroid_coordinates->data)[coord_index+1];
794      //compute x- and y- distances between the centroid of the FV triangle and that of its neighbour
795      dx1=x1-x; dy1=y1-y;
796      //set area2 as the square of the distance
797      area2=dx1*dx1+dy1*dy1;
798      //set dx2=(x1-x0)/((x1-x0)^2+(y1-y0)^2) and dy2=(y1-y0)/((x1-x0)^2+(y1-y0)^2) which
799      //respectively correspond to the x- and y- gradients of the conserved quantities
800      dx2=1.0/area2;
801      dy2=dx2*dy1;
802      dx2*=dx1;
803
804      //## stage ###
805      //compute differentials
806      dq1=((double *)stage_centroid_values->data)[k1]-((double *)stage_centroid_values->data)[k];
807      //calculate the gradient between the centroid of the FV triangle and that of its neighbour
808      a=dq1*dx2;
809      b=dq1*dy2;
810      //calculate provisional vertex jumps, to be limited
811      dqv[0]=a*dxv0+b*dyv0;
812      dqv[1]=a*dxv1+b*dyv1;
813      dqv[2]=a*dxv2+b*dyv2;
814      //now limit the jumps
815      if (dq1>=0.0){
816        qmin=0.0;
817        qmax=dq1;
818      }
819      else{
820        qmin=dq1;
821        qmax=0.0;
822      }
823      limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
824      for (i=0;i<3;i++)
825        ((double *)stage_vertex_values->data)[k3+i]=((double *)stage_centroid_values->data)[k]+dqv[i];
826
827      //## xmom ###
828      //compute differentials
829      dq1=((double *)xmom_centroid_values->data)[k1]-((double *)xmom_centroid_values->data)[k];
830      //calculate the gradient between the centroid of the FV triangle and that of its neighbour
831      a=dq1*dx2;
832      b=dq1*dy2;
833      //calculate provisional vertex jumps, to be limited
834      dqv[0]=a*dxv0+b*dyv0;
835      dqv[1]=a*dxv1+b*dyv1;
836      dqv[2]=a*dxv2+b*dyv2;
837      //now limit the jumps
838      if (dq1>=0.0){
839        qmin=0.0;
840        qmax=dq1;
841      }
842      else{
843        qmin=dq1;
844        qmax=0.0;
845      }
846      limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
847      for (i=0;i<3;i++)
848        ((double *)xmom_vertex_values->data)[k3+i]=((double *)xmom_centroid_values->data)[k]+dqv[i];
849
850      //## ymom ###
851      //compute differentials
852      dq1=((double *)ymom_centroid_values->data)[k1]-((double *)ymom_centroid_values->data)[k];
853      //calculate the gradient between the centroid of the FV triangle and that of its neighbour
854      a=dq1*dx2;
855      b=dq1*dy2;
856      //calculate provisional vertex jumps, to be limited
857      dqv[0]=a*dxv0+b*dyv0;
858      dqv[1]=a*dxv1+b*dyv1;
859      dqv[2]=a*dxv2+b*dyv2;
860      //now limit the jumps
861      if (dq1>=0.0){
862        qmin=0.0;
863        qmax=dq1;
864      }
865      else{
866        qmin=dq1;
867        qmax=0.0;
868      }
869      limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
870      for (i=0;i<3;i++)
871        ((double *)ymom_vertex_values->data)[k3+i]=((double *)ymom_centroid_values->data)[k]+dqv[i];
872    }//else [number_of_boudaries==2]
873  }//for k=0 to number_of_elements-1
874  return Py_BuildValue("");
875}//extrapolate_second-order_sw
876
877PyObject *rotate(PyObject *self, PyObject *args, PyObject *kwargs) {
878  //
879  // r = rotate(q, normal, direction=1)
880  //
881  // Where q is assumed to be a Float numeric array of length 3 and
882  // normal a Float numeric array of length 2.
883
884
885  PyObject *Q, *Normal;
886  PyArrayObject *q, *r, *normal;
887
888  static char *argnames[] = {"q", "normal", "direction", NULL};
889  int dimensions[1], i, direction=1;
890  double n1, n2;
891
892  // Convert Python arguments to C
893  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OO|i", argnames,
894                                   &Q, &Normal, &direction))
895    return NULL;
896
897  //Input checks (convert sequences into numeric arrays)
898  q = (PyArrayObject *)
899    PyArray_ContiguousFromObject(Q, PyArray_DOUBLE, 0, 0);
900  normal = (PyArrayObject *)
901    PyArray_ContiguousFromObject(Normal, PyArray_DOUBLE, 0, 0);
902
903
904  if (normal -> dimensions[0] != 2) {
905    PyErr_SetString(PyExc_RuntimeError, "Normal vector must have 2 components");
906    return NULL;
907  }
908
909  //Allocate space for return vector r (don't DECREF)
910  dimensions[0] = 3;
911  r = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
912
913  //Copy
914  for (i=0; i<3; i++) {
915    ((double *) (r -> data))[i] = ((double *) (q -> data))[i];
916  }
917
918  //Get normal and direction
919  n1 = ((double *) normal -> data)[0];
920  n2 = ((double *) normal -> data)[1];
921  if (direction == -1) n2 = -n2;
922
923  //Rotate
924  _rotate((double *) r -> data, n1, n2);
925
926  //Release numeric arrays
927  Py_DECREF(q);
928  Py_DECREF(normal);
929
930  //return result using PyArray to avoid memory leak
931  return PyArray_Return(r);
932}
933
934PyObject *compute_fluxes(PyObject *self, PyObject *args) {
935  /*Compute all fluxes and the timestep suitable for all volumes
936    in domain.
937
938    Compute total flux for each conserved quantity using "flux_function"
939
940    Fluxes across each edge are scaled by edgelengths and summed up
941    Resulting flux is then scaled by area and stored in
942    explicit_update for each of the three conserved quantities
943    stage, xmomentum and ymomentum
944
945    The maximal allowable speed computed by the flux_function for each volume
946    is converted to a timestep that must not be exceeded. The minimum of
947    those is computed as the next overall timestep.
948
949    Python call:
950    domain.timestep = compute_fluxes(timestep,
951                                     domain.epsilon,
952                                     domain.g,
953                                     domain.neighbours,
954                                     domain.neighbour_edges,
955                                     domain.normals,
956                                     domain.edgelengths,
957                                     domain.radii,
958                                     domain.areas,
959                                     Stage.edge_values,
960                                     Xmom.edge_values,
961                                     Ymom.edge_values,
962                                     Bed.edge_values,
963                                     Stage.boundary_values,
964                                     Xmom.boundary_values,
965                                     Ymom.boundary_values,
966                                     Stage.explicit_update,
967                                     Xmom.explicit_update,
968                                     Ymom.explicit_update,
969                                     already_computed_flux)
970
971
972    Post conditions:
973      domain.explicit_update is reset to computed flux values
974      domain.timestep is set to the largest step satisfying all volumes.
975
976
977  */
978
979
980  PyArrayObject *neighbours, *neighbour_edges,
981    *normals, *edgelengths, *radii, *areas,
982    *tri_full_flag,
983    *stage_edge_values,
984    *xmom_edge_values,
985    *ymom_edge_values,
986    *bed_edge_values,
987    *stage_boundary_values,
988    *xmom_boundary_values,
989    *ymom_boundary_values,
990    *stage_explicit_update,
991    *xmom_explicit_update,
992    *ymom_explicit_update,
993    *already_computed_flux;//tracks whether the flux across an edge has already been computed
994
995
996  //Local variables
997  double timestep, max_speed, epsilon, g;
998  double normal[2], ql[3], qr[3], zl, zr;
999  double edgeflux[3]; //Work arrays for summing up fluxes
1000
1001  int number_of_elements, k, i, m, n;
1002  int ki, nm=0, ki2; //Index shorthands
1003  static long call=1;
1004
1005
1006  // Convert Python arguments to C
1007  if (!PyArg_ParseTuple(args, "dddOOOOOOOOOOOOOOOOOO",
1008                        &timestep,
1009                        &epsilon,
1010                        &g,
1011                        &neighbours,
1012                        &neighbour_edges,
1013                        &normals,
1014                        &edgelengths, &radii, &areas,
1015                        &tri_full_flag,
1016                        &stage_edge_values,
1017                        &xmom_edge_values,
1018                        &ymom_edge_values,
1019                        &bed_edge_values,
1020                        &stage_boundary_values,
1021                        &xmom_boundary_values,
1022                        &ymom_boundary_values,
1023                        &stage_explicit_update,
1024                        &xmom_explicit_update,
1025                        &ymom_explicit_update,
1026                        &already_computed_flux)) {
1027    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
1028    return NULL;
1029  }
1030  number_of_elements = stage_edge_values -> dimensions[0];
1031  call++;//a static local variable to which already_computed_flux is compared
1032  //set explicit_update to zero for all conserved_quantities.
1033  //This assumes compute_fluxes called before forcing terms
1034  for (k=0; k<number_of_elements; k++) {
1035    ((double *) stage_explicit_update -> data)[k]=0.0;
1036    ((double *) xmom_explicit_update -> data)[k]=0.0;
1037    ((double *) ymom_explicit_update -> data)[k]=0.0;
1038  }
1039  //Loop through neighbours and compute edge flux for each
1040  for (k=0; k<number_of_elements; k++) {
1041    for (i=0; i<3; i++) {
1042      ki = k*3+i;
1043      if (((long *) already_computed_flux->data)[ki]==call)//we've already computed the flux across this edge
1044        continue;
1045      ql[0] = ((double *) stage_edge_values -> data)[ki];
1046      ql[1] = ((double *) xmom_edge_values -> data)[ki];
1047      ql[2] = ((double *) ymom_edge_values -> data)[ki];
1048      zl =    ((double *) bed_edge_values -> data)[ki];
1049
1050      //Quantities at neighbour on nearest face
1051      n = ((long *) neighbours -> data)[ki];
1052      if (n < 0) {
1053        m = -n-1; //Convert negative flag to index
1054        qr[0] = ((double *) stage_boundary_values -> data)[m];
1055        qr[1] = ((double *) xmom_boundary_values -> data)[m];
1056        qr[2] = ((double *) ymom_boundary_values -> data)[m];
1057        zr = zl; //Extend bed elevation to boundary
1058      } else {
1059        m = ((long *) neighbour_edges -> data)[ki];
1060        nm = n*3+m;
1061        qr[0] = ((double *) stage_edge_values -> data)[nm];
1062        qr[1] = ((double *) xmom_edge_values -> data)[nm];
1063        qr[2] = ((double *) ymom_edge_values -> data)[nm];
1064        zr =    ((double *) bed_edge_values -> data)[nm];
1065      }
1066      // Outward pointing normal vector
1067      // normal = domain.normals[k, 2*i:2*i+2]
1068      ki2 = 2*ki; //k*6 + i*2
1069      normal[0] = ((double *) normals -> data)[ki2];
1070      normal[1] = ((double *) normals -> data)[ki2+1];
1071      //Edge flux computation
1072      flux_function(ql, qr, zl, zr,
1073                    normal[0], normal[1],
1074                    epsilon, g,
1075                    edgeflux, &max_speed);
1076      //update triangle k
1077      ((long *) already_computed_flux->data)[ki]=call;
1078      ((double *) stage_explicit_update -> data)[k] -= edgeflux[0]*((double *) edgelengths -> data)[ki];
1079      ((double *) xmom_explicit_update -> data)[k] -= edgeflux[1]*((double *) edgelengths -> data)[ki];
1080      ((double *) ymom_explicit_update -> data)[k] -= edgeflux[2]*((double *) edgelengths -> data)[ki];
1081      //update the neighbour n
1082      if (n>=0){
1083        ((long *) already_computed_flux->data)[nm]=call;
1084        ((double *) stage_explicit_update -> data)[n] += edgeflux[0]*((double *) edgelengths -> data)[nm];
1085        ((double *) xmom_explicit_update -> data)[n] += edgeflux[1]*((double *) edgelengths -> data)[nm];
1086        ((double *) ymom_explicit_update -> data)[n] += edgeflux[2]*((double *) edgelengths -> data)[nm];
1087      }
1088      ///for (j=0; j<3; j++) {
1089        ///flux[j] -= edgeflux[j]*((double *) edgelengths -> data)[ki];
1090        ///}
1091        //Update timestep
1092        //timestep = min(timestep, domain.radii[k]/max_speed)
1093        //FIXME: SR Add parameter for CFL condition
1094    if ( ((long *) tri_full_flag -> data)[k] == 1) {
1095            if (max_speed > epsilon) {
1096                timestep = min(timestep, ((double *) radii -> data)[k]/max_speed);
1097                //maxspeed in flux_function is calculated as max(|u+a|,|u-a|)
1098                if (n>=0)
1099                    timestep = min(timestep, ((double *) radii -> data)[n]/max_speed);
1100            }
1101    }
1102    } // end for i
1103    //Normalise by area and store for when all conserved
1104    //quantities get updated
1105    ((double *) stage_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
1106    ((double *) xmom_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
1107    ((double *) ymom_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
1108  } //end for k
1109  return Py_BuildValue("d", timestep);
1110}
1111
1112PyObject *protect(PyObject *self, PyObject *args) {
1113  //
1114  //    protect(minimum_allowed_height, maximum_allowed_speed, wc, zc, xmomc, ymomc)
1115
1116
1117  PyArrayObject
1118  *wc,            //Stage at centroids
1119  *zc,            //Elevation at centroids
1120  *xmomc,         //Momentums at centroids
1121  *ymomc;
1122
1123
1124  int N;
1125  double minimum_allowed_height, maximum_allowed_speed, epsilon;
1126
1127  // Convert Python arguments to C
1128  if (!PyArg_ParseTuple(args, "dddOOOO",
1129                        &minimum_allowed_height,
1130                        &maximum_allowed_speed,
1131                        &epsilon,
1132                        &wc, &zc, &xmomc, &ymomc))
1133    return NULL;
1134
1135  N = wc -> dimensions[0];
1136
1137  _protect(N,
1138           minimum_allowed_height,
1139           maximum_allowed_speed,
1140           epsilon,
1141           (double*) wc -> data,
1142           (double*) zc -> data,
1143           (double*) xmomc -> data,
1144           (double*) ymomc -> data);
1145
1146  return Py_BuildValue("");
1147}
1148
1149
1150
1151PyObject *balance_deep_and_shallow(PyObject *self, PyObject *args) {
1152  //
1153  //    balance_deep_and_shallow(wc, zc, hc, wv, zv, hv,
1154  //                             xmomc, ymomc, xmomv, ymomv)
1155
1156
1157  PyArrayObject
1158    *wc,            //Stage at centroids
1159    *zc,            //Elevation at centroids
1160    *hc,            //Height at centroids
1161    *wv,            //Stage at vertices
1162    *zv,            //Elevation at vertices
1163    *hv,            //Depths at vertices
1164    *hvbar,         //h-Limited depths at vertices
1165    *xmomc,         //Momentums at centroids and vertices
1166    *ymomc,
1167    *xmomv,
1168    *ymomv;
1169
1170  int N; //, err;
1171
1172  // Convert Python arguments to C
1173  if (!PyArg_ParseTuple(args, "OOOOOOOOOOO",
1174                        &wc, &zc, &hc,
1175                        &wv, &zv, &hv, &hvbar,
1176                        &xmomc, &ymomc, &xmomv, &ymomv))
1177    return NULL;
1178
1179  N = wc -> dimensions[0];
1180
1181  _balance_deep_and_shallow(N,
1182                            (double*) wc -> data,
1183                            (double*) zc -> data,
1184                            (double*) hc -> data,
1185                            (double*) wv -> data,
1186                            (double*) zv -> data,
1187                            (double*) hv -> data,
1188                            (double*) hvbar -> data,
1189                            (double*) xmomc -> data,
1190                            (double*) ymomc -> data,
1191                            (double*) xmomv -> data,
1192                            (double*) ymomv -> data);
1193
1194
1195  return Py_BuildValue("");
1196}
1197
1198
1199
1200PyObject *h_limiter(PyObject *self, PyObject *args) {
1201
1202  PyObject *domain, *Tmp;
1203  PyArrayObject
1204    *hv, *hc, //Depth at vertices and centroids
1205    *hvbar,   //Limited depth at vertices (return values)
1206    *neighbours;
1207
1208  int k, i, n, N, k3;
1209  int dimensions[2];
1210  double beta_h; //Safety factor (see config.py)
1211  double *hmin, *hmax, hn;
1212
1213  // Convert Python arguments to C
1214  if (!PyArg_ParseTuple(args, "OOO", &domain, &hc, &hv))
1215    return NULL;
1216
1217  neighbours = get_consecutive_array(domain, "neighbours");
1218
1219  //Get safety factor beta_h
1220  Tmp = PyObject_GetAttrString(domain, "beta_h");
1221  if (!Tmp)
1222    return NULL;
1223
1224  beta_h = PyFloat_AsDouble(Tmp);
1225
1226  Py_DECREF(Tmp);
1227
1228  N = hc -> dimensions[0];
1229
1230  //Create hvbar
1231  dimensions[0] = N;
1232  dimensions[1] = 3;
1233  hvbar = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
1234
1235
1236  //Find min and max of this and neighbour's centroid values
1237  hmin = malloc(N * sizeof(double));
1238  hmax = malloc(N * sizeof(double));
1239  for (k=0; k<N; k++) {
1240    k3=k*3;
1241
1242    hmin[k] = ((double*) hc -> data)[k];
1243    hmax[k] = hmin[k];
1244
1245    for (i=0; i<3; i++) {
1246      n = ((long*) neighbours -> data)[k3+i];
1247
1248      //Initialise hvbar with values from hv
1249      ((double*) hvbar -> data)[k3+i] = ((double*) hv -> data)[k3+i];
1250
1251      if (n >= 0) {
1252        hn = ((double*) hc -> data)[n]; //Neighbour's centroid value
1253
1254        hmin[k] = min(hmin[k], hn);
1255        hmax[k] = max(hmax[k], hn);
1256      }
1257    }
1258  }
1259
1260  // Call underlying standard routine
1261  _limit(N, beta_h, (double*) hc -> data, (double*) hvbar -> data, hmin, hmax);
1262
1263  // // //Py_DECREF(domain); //FIXME: NEcessary?
1264  free(hmin);
1265  free(hmax);
1266
1267  //return result using PyArray to avoid memory leak
1268  return PyArray_Return(hvbar);
1269  //return Py_BuildValue("");
1270}
1271
1272PyObject *h_limiter_sw(PyObject *self, PyObject *args) {
1273  //a faster version of h_limiter above
1274  PyObject *domain, *Tmp;
1275  PyArrayObject
1276    *hv, *hc, //Depth at vertices and centroids
1277    *hvbar,   //Limited depth at vertices (return values)
1278    *neighbours;
1279
1280  int k, i, N, k3,k0,k1,k2;
1281  int dimensions[2];
1282  double beta_h; //Safety factor (see config.py)
1283  double hmin, hmax, dh[3];
1284  // Convert Python arguments to C
1285  if (!PyArg_ParseTuple(args, "OOO", &domain, &hc, &hv))
1286    return NULL;
1287  neighbours = get_consecutive_array(domain, "neighbours");
1288
1289  //Get safety factor beta_h
1290  Tmp = PyObject_GetAttrString(domain, "beta_h");
1291  if (!Tmp)
1292    return NULL;
1293  beta_h = PyFloat_AsDouble(Tmp);
1294
1295  Py_DECREF(Tmp);
1296
1297  N = hc -> dimensions[0];
1298
1299  //Create hvbar
1300  dimensions[0] = N;
1301  dimensions[1] = 3;
1302  hvbar = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
1303  for (k=0;k<N;k++){
1304    k3=k*3;
1305    //get the ids of the neighbours
1306    k0 = ((long*) neighbours -> data)[k3];
1307    k1 = ((long*) neighbours -> data)[k3+1];
1308    k2 = ((long*) neighbours -> data)[k3+2];
1309    //set hvbar provisionally
1310    for (i=0;i<3;i++){
1311      ((double*) hvbar -> data)[k3+i] = ((double*) hv -> data)[k3+i];
1312      dh[i]=((double*) hvbar -> data)[k3+i]-((double*) hc -> data)[k];
1313    }
1314    hmin=((double*) hc -> data)[k];
1315    hmax=hmin;
1316    if (k0>=0){
1317      hmin=min(hmin,((double*) hc -> data)[k0]);
1318      hmax=max(hmax,((double*) hc -> data)[k0]);
1319    }
1320    if (k1>=0){
1321      hmin=min(hmin,((double*) hc -> data)[k1]);
1322      hmax=max(hmax,((double*) hc -> data)[k1]);
1323    }
1324    if (k2>=0){
1325      hmin=min(hmin,((double*) hc -> data)[k2]);
1326      hmax=max(hmax,((double*) hc -> data)[k2]);
1327    }
1328    hmin-=((double*) hc -> data)[k];
1329    hmax-=((double*) hc -> data)[k];
1330    limit_gradient(dh,hmin,hmax,beta_h);
1331    for (i=0;i<3;i++)
1332      ((double*) hvbar -> data)[k3+i] = ((double*) hc -> data)[k]+dh[i];
1333  }
1334  return PyArray_Return(hvbar);
1335}
1336
1337PyObject *assign_windfield_values(PyObject *self, PyObject *args) {
1338  //
1339  //      assign_windfield_values(xmom_update, ymom_update,
1340  //                              s_vec, phi_vec, self.const)
1341
1342
1343
1344  PyArrayObject   //(one element per triangle)
1345  *s_vec,         //Speeds
1346  *phi_vec,       //Bearings
1347  *xmom_update,   //Momentum updates
1348  *ymom_update;
1349
1350
1351  int N;
1352  double cw;
1353
1354  // Convert Python arguments to C
1355  if (!PyArg_ParseTuple(args, "OOOOd",
1356                        &xmom_update,
1357                        &ymom_update,
1358                        &s_vec, &phi_vec,
1359                        &cw))
1360    return NULL;
1361
1362  N = xmom_update -> dimensions[0];
1363
1364  _assign_wind_field_values(N,
1365           (double*) xmom_update -> data,
1366           (double*) ymom_update -> data,
1367           (double*) s_vec -> data,
1368           (double*) phi_vec -> data,
1369           cw);
1370
1371  return Py_BuildValue("");
1372}
1373
1374
1375
1376
1377//////////////////////////////////////////
1378// Method table for python module
1379static struct PyMethodDef MethodTable[] = {
1380  /* The cast of the function is necessary since PyCFunction values
1381   * only take two PyObject* parameters, and rotate() takes
1382   * three.
1383   */
1384
1385  {"rotate", (PyCFunction)rotate, METH_VARARGS | METH_KEYWORDS, "Print out"},
1386  {"extrapolate_second_order_sw", extrapolate_second_order_sw, METH_VARARGS, "Print out"},
1387  {"compute_fluxes", compute_fluxes, METH_VARARGS, "Print out"},
1388  {"gravity", gravity, METH_VARARGS, "Print out"},
1389  {"manning_friction", manning_friction, METH_VARARGS, "Print out"},
1390  {"balance_deep_and_shallow", balance_deep_and_shallow,
1391   METH_VARARGS, "Print out"},
1392  {"h_limiter", h_limiter,
1393   METH_VARARGS, "Print out"},
1394  {"h_limiter_sw", h_limiter_sw,
1395   METH_VARARGS, "Print out"},
1396  {"protect", protect, METH_VARARGS | METH_KEYWORDS, "Print out"},
1397  {"assign_windfield_values", assign_windfield_values,
1398   METH_VARARGS | METH_KEYWORDS, "Print out"},
1399  //{"distribute_to_vertices_and_edges",
1400  // distribute_to_vertices_and_edges, METH_VARARGS},
1401  //{"update_conserved_quantities",
1402  // update_conserved_quantities, METH_VARARGS},
1403  //{"set_initialcondition",
1404  // set_initialcondition, METH_VARARGS},
1405  {NULL, NULL}
1406};
1407
1408// Module initialisation
1409void initshallow_water_ext(void){
1410  Py_InitModule("shallow_water_ext", MethodTable);
1411
1412  import_array();     //Necessary for handling of NumPY structures
1413}
Note: See TracBrowser for help on using the repository browser.