source: inundation/pyvolution/shallow_water_ext.c @ 2567

Last change on this file since 2567 was 2567, checked in by ole, 18 years ago

Cleanup and comments

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