source:inundation/pyvolution/shallow_water_ext.c@2621

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

File size: 42.0 KB
Line
1// Python - C extension module for shallow_water.py
2//
3// To compile (Python2.3):
4//  gcc -c domain_ext.c -I/usr/include/python2.3 -o domain_ext.o -Wall -O
5//  gcc -shared domain_ext.o  -o domain_ext.so
6//
7// or use python compile.py
8//
9// See the module shallow_water.py
10//
11//
12// Ole Nielsen, GA 2004
13
14
15#include "Python.h"
16#include "Numeric/arrayobject.h"
17#include "math.h"
18#include <stdio.h>
19
20//Shared code snippets
21#include "util_ext.h"
22
23const double pi = 3.14159265358979;
24
25// Computational function for rotation
26int _rotate(double *q, double n1, double n2) {
27  /*Rotate the momentum component q (q[1], q[2])
28    from x,y coordinates to coordinates based on normal vector (n1, n2).
29
30    Result is returned in array 3x1 r
31    To rotate in opposite direction, call rotate with (q, n1, -n2)
32
33    Contents of q are changed by this function */
34
35
36  double q1, q2;
37
38  //Shorthands
39  q1 = q[1];  //uh momentum
40  q2 = q[2];  //vh momentum
41
42  //Rotate
43  q[1] =  n1*q1 + n2*q2;
44  q[2] = -n2*q1 + n1*q2;
45
46  return 0;
47}
48
49int find_qmin_and_qmax(double dq0, double dq1, double dq2, double *qmin, double *qmax){
50  //Considering the centroid of an FV triangle and the vertices of its auxiliary triangle, find
51  //qmin=min(q)-qc and qmax=max(q)-qc, where min(q) and max(q) are respectively min and max over the
52  //four values (at the centroid of the FV triangle and the auxiliary triangle vertices),
53  //and qc is the centroid
54  //dq0=q(vertex0)-q(centroid of FV triangle)
55  //dq1=q(vertex1)-q(vertex0)
56  //dq2=q(vertex2)-q(vertex0)
57  if (dq0>=0.0){
58    if (dq1>=dq2){
59      if (dq1>=0.0)
60        *qmax=dq0+dq1;
61      else
62        *qmax=dq0;
63      if ((*qmin=dq0+dq2)<0)
64        ;//qmin is already set to correct value
65      else
66        *qmin=0.0;
67    }
68    else{//dq1<dq2
69      if (dq2>0)
70        *qmax=dq0+dq2;
71      else
72        *qmax=dq0;
73      if ((*qmin=dq0+dq1)<0)
74        ;//qmin is the correct value
75      else
76        *qmin=0.0;
77    }
78  }
79  else{//dq0<0
80    if (dq1<=dq2){
81      if (dq1<0.0)
82        *qmin=dq0+dq1;
83      else
84        *qmin=dq0;
85      if ((*qmax=dq0+dq2)>0.0)
86        ;//qmax is already set to the correct value
87      else
88        *qmax=0.0;
89    }
90    else{//dq1>dq2
91      if (dq2<0.0)
92        *qmin=dq0+dq2;
93      else
94        *qmin=dq0;
95      if ((*qmax=dq0+dq1)>0.0)
96        ;//qmax is already set to the correct value
97      else
98        *qmax=0.0;
99    }
100  }
101  return 0;
102}
103
104int limit_gradient(double *dqv, double qmin, double qmax, double beta_w){
105  //given provisional jumps dqv from the FV triangle centroid to its vertices and
106  //jumps qmin (qmax) between the centroid of the FV triangle and the
107  //minimum (maximum) of the values at the centroid of the FV triangle and the auxiliary triangle vertices,
108  //calculate a multiplicative factor phi by which the provisional vertex jumps are to be limited
109  int i;
110  double r=1000.0, r0=1.0, phi=1.0;
111  static double TINY = 1.0e-100;//to avoid machine accuracy problems.
112  //Any provisional jump with magnitude < TINY does not contribute to the limiting process.
113  for (i=0;i<3;i++){
114    if (dqv[i]<-TINY)
115      r0=qmin/dqv[i];
116    if (dqv[i]>TINY)
117      r0=qmax/dqv[i];
118    r=min(r0,r);
119    //
120  }
121  phi=min(r*beta_w,1.0);
122  for (i=0;i<3;i++)
123    dqv[i]=dqv[i]*phi;
124  return 0;
125}
126
127// 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      //Old code
399      //wc[k] = zc[k]; //Contain 'lost mass' error
400      //printf("Speed has been reduced to %.3f\n", 0.0);
401      //xmomc[k] = 0.0;
402      //ymomc[k] = 0.0;
403    }
404  }
405  return 0;
406}
407
408
409
410
411int _assign_wind_field_values(int N,
412                              double* xmom_update,
413                              double* ymom_update,
414                              double* s_vec,
415                              double* phi_vec,
416                              double cw) {
417
418  //Assign windfield values to momentum updates
419
420  int k;
421  double S, s, phi, u, v;
422
423  for (k=0; k<N; k++) {
424
425    s = s_vec[k];
426    phi = phi_vec[k];
427
429    phi = phi*pi/180;
430
431    //Compute velocity vector (u, v)
432    u = s*cos(phi);
433    v = s*sin(phi);
434
435    //Compute wind stress
436    S = cw * sqrt(u*u + v*v);
437    xmom_update[k] += S*u;
438    ymom_update[k] += S*v;
439  }
440  return 0;
441}
442
443
444
445///////////////////////////////////////////////////////////////////
446// Gateways to Python
447
448PyObject *gravity(PyObject *self, PyObject *args) {
449  //
450  //  gravity(g, h, v, x, xmom, ymom)
451  //
452
453
454  PyArrayObject *h, *v, *x, *xmom, *ymom;
455  int k, i, N, k3, k6;
456  double g, avg_h, zx, zy;
457  double x0, y0, x1, y1, x2, y2, z0, z1, z2;
458
459  if (!PyArg_ParseTuple(args, "dOOOOO",
460                        &g, &h, &v, &x,
461                        &xmom, &ymom))
462    return NULL;
463
464  N = h -> dimensions[0];
465  for (k=0; k<N; k++) {
466    k3 = 3*k;  // base index
467    k6 = 6*k;  // base index
468
469    avg_h = 0.0;
470    for (i=0; i<3; i++) {
471      avg_h += ((double *) h -> data)[k3+i];
472    }
473    avg_h /= 3;
474
475
476    //Compute bed slope
477    x0 = ((double*) x -> data)[k6 + 0];
478    y0 = ((double*) x -> data)[k6 + 1];
479    x1 = ((double*) x -> data)[k6 + 2];
480    y1 = ((double*) x -> data)[k6 + 3];
481    x2 = ((double*) x -> data)[k6 + 4];
482    y2 = ((double*) x -> data)[k6 + 5];
483
484
485    z0 = ((double*) v -> data)[k3 + 0];
486    z1 = ((double*) v -> data)[k3 + 1];
487    z2 = ((double*) v -> data)[k3 + 2];
488
489    _gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2, &zx, &zy);
490
491    //Update momentum
492    ((double*) xmom -> data)[k] += -g*zx*avg_h;
493    ((double*) ymom -> data)[k] += -g*zy*avg_h;
494  }
495
496  return Py_BuildValue("");
497}
498
499
500PyObject *manning_friction(PyObject *self, PyObject *args) {
501  //
502  // manning_friction(g, eps, h, uh, vh, eta, xmom_update, ymom_update)
503  //
504
505
506  PyArrayObject *w, *z, *uh, *vh, *eta, *xmom, *ymom;
507  int N;
508  double g, eps;
509
510  if (!PyArg_ParseTuple(args, "ddOOOOOOO",
511                        &g, &eps, &w, &z, &uh, &vh, &eta,
512                        &xmom, &ymom))
513    return NULL;
514
515  N = w -> dimensions[0];
516  _manning_friction(g, eps, N,
517                    (double*) w -> data,
518                    (double*) z -> data,
519                    (double*) uh -> data,
520                    (double*) vh -> data,
521                    (double*) eta -> data,
522                    (double*) xmom -> data,
523                    (double*) ymom -> data);
524
525  return Py_BuildValue("");
526}
527
528PyObject *extrapolate_second_order_sw(PyObject *self, PyObject *args) {
529  /*Compute the vertex values based on a linear reconstruction on each triangle
530    These values are calculated as follows:
531    1) For each triangle not adjacent to a boundary, we consider the auxiliary triangle
532    formed by the centroids of its three neighbours.
533    2) For each conserved quantity, we integrate around the auxiliary triangle's boundary the product
534    of the quantity and the outward normal vector. Dividing by the triangle area gives (a,b), the average
535    of the vector (q_x,q_y) on the auxiliary triangle. We suppose that the linear reconstruction on the
536    original triangle has gradient (a,b).
537    3) Provisional vertex junmps dqv[0,1,2] are computed and these are then limited by calling the functions
539
540    Python call:
541    extrapolate_second_order_sw(domain.surrogate_neighbours,
542                                domain.number_of_boundaries
543                                domain.centroid_coordinates,
544                                Stage.centroid_values
545                                Xmom.centroid_values
546                                Ymom.centroid_values
547                                domain.vertex_coordinates,
548                                Stage.vertex_values,
549                                Xmom.vertex_values,
550                                Ymom.vertex_values)
551
552    Post conditions:
553            The vertices of each triangle have values from a limited linear reconstruction
554            based on centroid values
555
556  */
557  PyArrayObject *surrogate_neighbours,
558    *number_of_boundaries,
559    *centroid_coordinates,
560    *stage_centroid_values,
561    *xmom_centroid_values,
562    *ymom_centroid_values,
563    *vertex_coordinates,
564    *stage_vertex_values,
565    *xmom_vertex_values,
566    *ymom_vertex_values;
567  PyObject *domain, *Tmp;
568  //Local variables
569  double a, b;//gradient vector, not stored but used to calculate vertex values from centroids
570  int number_of_elements,k,k0,k1,k2,k3,k6,coord_index,i;
571  double x,y,x0,y0,x1,y1,x2,y2,xv0,yv0,xv1,yv1,xv2,yv2;//vertices of the auxiliary triangle
572  double dx1,dx2,dy1,dy2,dxv0,dxv1,dxv2,dyv0,dyv1,dyv2,dq0,dq1,dq2,area2;
573  double dqv[3], qmin, qmax, beta_w;//provisional jumps from centroids to v'tices and safety factor re limiting
574  //by which these jumps are limited
575  // Convert Python arguments to C
576  if (!PyArg_ParseTuple(args, "OOOOOOOOOOO",
577                        &domain,
578                        &surrogate_neighbours,
579                        &number_of_boundaries,
580                        &centroid_coordinates,
581                        &stage_centroid_values,
582                        &xmom_centroid_values,
583                        &ymom_centroid_values,
584                        &vertex_coordinates,
585                        &stage_vertex_values,
586                        &xmom_vertex_values,
587                        &ymom_vertex_values)) {
588    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
589    return NULL;
590  }
591
592  //get the safety factor beta_w, set in the config.py file. This is used in the limiting process
593  Tmp = PyObject_GetAttrString(domain, "beta_w");
594  if (!Tmp)
595    return NULL;
596  beta_w = PyFloat_AsDouble(Tmp);
597  Py_DECREF(Tmp);
598  number_of_elements = stage_centroid_values -> dimensions[0];
599  for (k=0; k<number_of_elements; k++) {
600    k3=k*3;
601    k6=k*6;
602
603    if (((long *) number_of_boundaries->data)[k]==3){/*no neighbours, set gradient on the triangle to zero*/
604      ((double *) stage_vertex_values->data)[k3]=((double *)stage_centroid_values->data)[k];
605      ((double *) stage_vertex_values->data)[k3+1]=((double *)stage_centroid_values->data)[k];
606      ((double *) stage_vertex_values->data)[k3+2]=((double *)stage_centroid_values->data)[k];
607      ((double *) xmom_vertex_values->data)[k3]=((double *)xmom_centroid_values->data)[k];
608      ((double *) xmom_vertex_values->data)[k3+1]=((double *)xmom_centroid_values->data)[k];
609      ((double *) xmom_vertex_values->data)[k3+2]=((double *)xmom_centroid_values->data)[k];
610      ((double *) ymom_vertex_values->data)[k3]=((double *)ymom_centroid_values->data)[k];
611      ((double *) ymom_vertex_values->data)[k3+1]=((double *)ymom_centroid_values->data)[k];
612      ((double *) ymom_vertex_values->data)[k3+2]=((double *)ymom_centroid_values->data)[k];
613      continue;
614    }
615    else {//we will need centroid coordinates and vertex coordinates of the triangle
616      //get the vertex coordinates of the FV triangle
617      xv0=((double *)vertex_coordinates->data)[k6]; yv0=((double *)vertex_coordinates->data)[k6+1];
618      xv1=((double *)vertex_coordinates->data)[k6+2]; yv1=((double *)vertex_coordinates->data)[k6+3];
619      xv2=((double *)vertex_coordinates->data)[k6+4]; yv2=((double *)vertex_coordinates->data)[k6+5];
620      //get the centroid coordinates of the FV triangle
621      coord_index=2*k;
622      x=((double *)centroid_coordinates->data)[coord_index];
623      y=((double *)centroid_coordinates->data)[coord_index+1];
624      //store x- and y- differentials for the vertices of the FV triangle relative to the centroid
625      dxv0=xv0-x; dxv1=xv1-x; dxv2=xv2-x;
626      dyv0=yv0-y; dyv1=yv1-y; dyv2=yv2-y;
627    }
628    if (((long *)number_of_boundaries->data)[k]<=1){
629      //if no boundaries, auxiliary triangle is formed from the centroids of the three neighbours
630      //if one boundary, auxiliary triangle is formed from this centroid and its two neighbours
631      k0=((long *)surrogate_neighbours->data)[k3];
632      k1=((long *)surrogate_neighbours->data)[k3+1];
633      k2=((long *)surrogate_neighbours->data)[k3+2];
634      //get the auxiliary triangle's vertex coordinates (really the centroids of neighbouring triangles)
635      coord_index=2*k0;
636      x0=((double *)centroid_coordinates->data)[coord_index];
637      y0=((double *)centroid_coordinates->data)[coord_index+1];
638      coord_index=2*k1;
639      x1=((double *)centroid_coordinates->data)[coord_index];
640      y1=((double *)centroid_coordinates->data)[coord_index+1];
641      coord_index=2*k2;
642      x2=((double *)centroid_coordinates->data)[coord_index];
643      y2=((double *)centroid_coordinates->data)[coord_index+1];
644      //store x- and y- differentials for the vertices of the auxiliary triangle
645      dx1=x1-x0; dx2=x2-x0;
646      dy1=y1-y0; dy2=y2-y0;
647      //calculate 2*area of the auxiliary triangle
648      area2 = dy2*dx1 - dy1*dx2;//the triangle is guaranteed to be counter-clockwise
649      //If the mesh is 'weird' near the boundary, the trianlge might be flat or clockwise:
650      if (area2<=0)
651        return NULL;
652
653      //### stage ###
654      //calculate the difference between vertex 0 of the auxiliary triangle and the FV triangle centroid
655      dq0=((double *)stage_centroid_values->data)[k0]-((double *)stage_centroid_values->data)[k];
656      //calculate differentials between the vertices of the auxiliary triangle
657      dq1=((double *)stage_centroid_values->data)[k1]-((double *)stage_centroid_values->data)[k0];
658      dq2=((double *)stage_centroid_values->data)[k2]-((double *)stage_centroid_values->data)[k0];
659      //calculate the gradient of stage on the auxiliary triangle
660      a = dy2*dq1 - dy1*dq2;
661      a /= area2;
662      b = dx1*dq2 - dx2*dq1;
663      b /= area2;
664      //calculate provisional jumps in stage from the centroid of the FV tri to its vertices, to be limited
665      dqv[0]=a*dxv0+b*dyv0;
666      dqv[1]=a*dxv1+b*dyv1;
667      dqv[2]=a*dxv2+b*dyv2;
668      //now we want to find min and max of the centroid and the vertices of the auxiliary triangle
669      //and compute jumps from the centroid to the min and max
670      find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
672      for (i=0;i<3;i++)
673        ((double *)stage_vertex_values->data)[k3+i]=((double *)stage_centroid_values->data)[k]+dqv[i];
674
675      //### xmom ###
676      //calculate the difference between vertex 0 of the auxiliary triangle and the FV triangle centroid
677      dq0=((double *)xmom_centroid_values->data)[k0]-((double *)xmom_centroid_values->data)[k];
678      //calculate differentials between the vertices of the auxiliary triangle
679      dq1=((double *)xmom_centroid_values->data)[k1]-((double *)xmom_centroid_values->data)[k0];
680      dq2=((double *)xmom_centroid_values->data)[k2]-((double *)xmom_centroid_values->data)[k0];
681      //calculate the gradient of xmom on the auxiliary triangle
682      a = dy2*dq1 - dy1*dq2;
683      a /= area2;
684      b = dx1*dq2 - dx2*dq1;
685      b /= area2;
686      //calculate provisional jumps in stage from the centroid of the FV tri to its vertices, to be limited
687      dqv[0]=a*dxv0+b*dyv0;
688      dqv[1]=a*dxv1+b*dyv1;
689      dqv[2]=a*dxv2+b*dyv2;
690      //now we want to find min and max of the centroid and the vertices of the auxiliary triangle
691      //and compute jumps from the centroid to the min and max
692      find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
694      for (i=0;i<3;i++)
695        ((double *)xmom_vertex_values->data)[k3+i]=((double *)xmom_centroid_values->data)[k]+dqv[i];
696
697      //### ymom ###
698      //calculate the difference between vertex 0 of the auxiliary triangle and the FV triangle centroid
699      dq0=((double *)ymom_centroid_values->data)[k0]-((double *)ymom_centroid_values->data)[k];
700      //calculate differentials between the vertices of the auxiliary triangle
701      dq1=((double *)ymom_centroid_values->data)[k1]-((double *)ymom_centroid_values->data)[k0];
702      dq2=((double *)ymom_centroid_values->data)[k2]-((double *)ymom_centroid_values->data)[k0];
703      //calculate the gradient of xmom on the auxiliary triangle
704      a = dy2*dq1 - dy1*dq2;
705      a /= area2;
706      b = dx1*dq2 - dx2*dq1;
707      b /= area2;
708      //calculate provisional jumps in stage from the centroid of the FV tri to its vertices, to be limited
709      dqv[0]=a*dxv0+b*dyv0;
710      dqv[1]=a*dxv1+b*dyv1;
711      dqv[2]=a*dxv2+b*dyv2;
712      //now we want to find min and max of the centroid and the vertices of the auxiliary triangle
713      //and compute jumps from the centroid to the min and max
714      find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
716      for (i=0;i<3;i++)
717        ((double *)ymom_vertex_values->data)[k3+i]=((double *)ymom_centroid_values->data)[k]+dqv[i];
718    }//if (number_of_boundaries[k]<=1)
719    else{//number_of_boundaries==2
720      //one internal neighbour and gradient is in direction of the neighbour's centroid
721      //find the only internal neighbour
722      for (k2=k3;k2<k3+3;k2++){//k2 just indexes the edges of triangle k
723        if (((long *)surrogate_neighbours->data)[k2]!=k)//find internal neighbour of triabngle k
724          break;
725      }
726      if ((k2==k3+3))//if we didn't find an internal neighbour
727        return NULL;//error
728      k1=((long *)surrogate_neighbours->data)[k2];
729      //the coordinates of the triangle are already (x,y). Get centroid of the neighbour (x1,y1)
730      coord_index=2*k1;
731      x1=((double *)centroid_coordinates->data)[coord_index];
732      y1=((double *)centroid_coordinates->data)[coord_index+1];
733      //compute x- and y- distances between the centroid of the FV triangle and that of its neighbour
734      dx1=x1-x; dy1=y1-y;
735      //set area2 as the square of the distance
736      area2=dx1*dx1+dy1*dy1;
737      //set dx2=(x1-x0)/((x1-x0)^2+(y1-y0)^2) and dy2=(y1-y0)/((x1-x0)^2+(y1-y0)^2) which
738      //respectively correspond to the x- and y- gradients of the conserved quantities
739      dx2=1.0/area2;
740      dy2=dx2*dy1;
741      dx2*=dx1;
742
743      //## stage ###
744      //compute differentials
745      dq1=((double *)stage_centroid_values->data)[k1]-((double *)stage_centroid_values->data)[k];
746      //calculate the gradient between the centroid of the FV triangle and that of its neighbour
747      a=dq1*dx2;
748      b=dq1*dy2;
749      //calculate provisional vertex jumps, to be limited
750      dqv[0]=a*dxv0+b*dyv0;
751      dqv[1]=a*dxv1+b*dyv1;
752      dqv[2]=a*dxv2+b*dyv2;
753      //now limit the jumps
754      if (dq1>=0.0){
755        qmin=0.0;
756        qmax=dq1;
757      }
758      else{
759        qmin=dq1;
760        qmax=0.0;
761      }
763      for (i=0;i<3;i++)
764        ((double *)stage_vertex_values->data)[k3+i]=((double *)stage_centroid_values->data)[k]+dqv[i];
765
766      //## xmom ###
767      //compute differentials
768      dq1=((double *)xmom_centroid_values->data)[k1]-((double *)xmom_centroid_values->data)[k];
769      //calculate the gradient between the centroid of the FV triangle and that of its neighbour
770      a=dq1*dx2;
771      b=dq1*dy2;
772      //calculate provisional vertex jumps, to be limited
773      dqv[0]=a*dxv0+b*dyv0;
774      dqv[1]=a*dxv1+b*dyv1;
775      dqv[2]=a*dxv2+b*dyv2;
776      //now limit the jumps
777      if (dq1>=0.0){
778        qmin=0.0;
779        qmax=dq1;
780      }
781      else{
782        qmin=dq1;
783        qmax=0.0;
784      }
786      for (i=0;i<3;i++)
787        ((double *)xmom_vertex_values->data)[k3+i]=((double *)xmom_centroid_values->data)[k]+dqv[i];
788
789      //## ymom ###
790      //compute differentials
791      dq1=((double *)ymom_centroid_values->data)[k1]-((double *)ymom_centroid_values->data)[k];
792      //calculate the gradient between the centroid of the FV triangle and that of its neighbour
793      a=dq1*dx2;
794      b=dq1*dy2;
795      //calculate provisional vertex jumps, to be limited
796      dqv[0]=a*dxv0+b*dyv0;
797      dqv[1]=a*dxv1+b*dyv1;
798      dqv[2]=a*dxv2+b*dyv2;
799      //now limit the jumps
800      if (dq1>=0.0){
801        qmin=0.0;
802        qmax=dq1;
803      }
804      else{
805        qmin=dq1;
806        qmax=0.0;
807      }
809      for (i=0;i<3;i++)
810        ((double *)ymom_vertex_values->data)[k3+i]=((double *)ymom_centroid_values->data)[k]+dqv[i];
811    }//else [number_of_boudaries==2]
812  }//for k=0 to number_of_elements-1
813  return Py_BuildValue("");
814}//extrapolate_second-order_sw
815
816PyObject *rotate(PyObject *self, PyObject *args, PyObject *kwargs) {
817  //
818  // r = rotate(q, normal, direction=1)
819  //
820  // Where q is assumed to be a Float numeric array of length 3 and
821  // normal a Float numeric array of length 2.
822
823
824  PyObject *Q, *Normal;
825  PyArrayObject *q, *r, *normal;
826
827  static char *argnames[] = {"q", "normal", "direction", NULL};
828  int dimensions[1], i, direction=1;
829  double n1, n2;
830
831  // Convert Python arguments to C
832  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OO|i", argnames,
833                                   &Q, &Normal, &direction))
834    return NULL;
835
836  //Input checks (convert sequences into numeric arrays)
837  q = (PyArrayObject *)
838    PyArray_ContiguousFromObject(Q, PyArray_DOUBLE, 0, 0);
839  normal = (PyArrayObject *)
840    PyArray_ContiguousFromObject(Normal, PyArray_DOUBLE, 0, 0);
841
842
843  if (normal -> dimensions[0] != 2) {
844    PyErr_SetString(PyExc_RuntimeError, "Normal vector must have 2 components");
845    return NULL;
846  }
847
848  //Allocate space for return vector r (don't DECREF)
849  dimensions[0] = 3;
850  r = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
851
852  //Copy
853  for (i=0; i<3; i++) {
854    ((double *) (r -> data))[i] = ((double *) (q -> data))[i];
855  }
856
857  //Get normal and direction
858  n1 = ((double *) normal -> data)[0];
859  n2 = ((double *) normal -> data)[1];
860  if (direction == -1) n2 = -n2;
861
862  //Rotate
863  _rotate((double *) r -> data, n1, n2);
864
865  //Release numeric arrays
866  Py_DECREF(q);
867  Py_DECREF(normal);
868
869  //return result using PyArray to avoid memory leak
870  return PyArray_Return(r);
871}
872
873PyObject *compute_fluxes(PyObject *self, PyObject *args) {
874  /*Compute all fluxes and the timestep suitable for all volumes
875    in domain.
876
877    Compute total flux for each conserved quantity using "flux_function"
878
879    Fluxes across each edge are scaled by edgelengths and summed up
880    Resulting flux is then scaled by area and stored in
881    explicit_update for each of the three conserved quantities
882    stage, xmomentum and ymomentum
883
884    The maximal allowable speed computed by the flux_function for each volume
885    is converted to a timestep that must not be exceeded. The minimum of
886    those is computed as the next overall timestep.
887
888    Python call:
889    domain.timestep = compute_fluxes(timestep,
890                                     domain.epsilon,
891                                     domain.g,
892                                     domain.neighbours,
893                                     domain.neighbour_edges,
894                                     domain.normals,
895                                     domain.edgelengths,
897                                     domain.areas,
898                                     Stage.edge_values,
899                                     Xmom.edge_values,
900                                     Ymom.edge_values,
901                                     Bed.edge_values,
902                                     Stage.boundary_values,
903                                     Xmom.boundary_values,
904                                     Ymom.boundary_values,
905                                     Stage.explicit_update,
906                                     Xmom.explicit_update,
907                                     Ymom.explicit_update,
909
910
911    Post conditions:
912      domain.explicit_update is reset to computed flux values
913      domain.timestep is set to the largest step satisfying all volumes.
914
915
916  */
917
918
919  PyArrayObject *neighbours, *neighbour_edges,
921    *stage_edge_values,
922    *xmom_edge_values,
923    *ymom_edge_values,
924    *bed_edge_values,
925    *stage_boundary_values,
926    *xmom_boundary_values,
927    *ymom_boundary_values,
928    *stage_explicit_update,
929    *xmom_explicit_update,
930    *ymom_explicit_update,
931    *already_computed_flux;//tracks whether the flux across an edge has already been computed
932
933
934  //Local variables
935  double timestep, max_speed, epsilon, g;
936  double normal[2], ql[3], qr[3], zl, zr;
937  double edgeflux[3]; //Work arrays for summing up fluxes
938
939  int number_of_elements, k, i, m, n;
940  int ki, nm=0, ki2; //Index shorthands
941  static long call=1;
942
943
944  // Convert Python arguments to C
945  if (!PyArg_ParseTuple(args, "dddOOOOOOOOOOOOOOOOO",
946                        &timestep,
947                        &epsilon,
948                        &g,
949                        &neighbours,
950                        &neighbour_edges,
951                        &normals,
953                        &stage_edge_values,
954                        &xmom_edge_values,
955                        &ymom_edge_values,
956                        &bed_edge_values,
957                        &stage_boundary_values,
958                        &xmom_boundary_values,
959                        &ymom_boundary_values,
960                        &stage_explicit_update,
961                        &xmom_explicit_update,
962                        &ymom_explicit_update,
964    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
965    return NULL;
966  }
967  number_of_elements = stage_edge_values -> dimensions[0];
968  call++;//a static local variable to which already_computed_flux is compared
969  //set explicit_update to zero for all conserved_quantities.
970  //This assumes compute_fluxes called before forcing terms
971  for (k=0; k<number_of_elements; k++) {
972    ((double *) stage_explicit_update -> data)[k]=0.0;
973    ((double *) xmom_explicit_update -> data)[k]=0.0;
974    ((double *) ymom_explicit_update -> data)[k]=0.0;
975  }
976  //Loop through neighbours and compute edge flux for each
977  for (k=0; k<number_of_elements; k++) {
978    for (i=0; i<3; i++) {
979      ki = k*3+i;
980      if (((long *) already_computed_flux->data)[ki]==call)//we've already computed the flux across this edge
981        continue;
982      ql[0] = ((double *) stage_edge_values -> data)[ki];
983      ql[1] = ((double *) xmom_edge_values -> data)[ki];
984      ql[2] = ((double *) ymom_edge_values -> data)[ki];
985      zl =    ((double *) bed_edge_values -> data)[ki];
986
987      //Quantities at neighbour on nearest face
988      n = ((long *) neighbours -> data)[ki];
989      if (n < 0) {
990        m = -n-1; //Convert negative flag to index
991        qr[0] = ((double *) stage_boundary_values -> data)[m];
992        qr[1] = ((double *) xmom_boundary_values -> data)[m];
993        qr[2] = ((double *) ymom_boundary_values -> data)[m];
994        zr = zl; //Extend bed elevation to boundary
995      } else {
996        m = ((long *) neighbour_edges -> data)[ki];
997        nm = n*3+m;
998        qr[0] = ((double *) stage_edge_values -> data)[nm];
999        qr[1] = ((double *) xmom_edge_values -> data)[nm];
1000        qr[2] = ((double *) ymom_edge_values -> data)[nm];
1001        zr =    ((double *) bed_edge_values -> data)[nm];
1002      }
1003      // Outward pointing normal vector
1004      // normal = domain.normals[k, 2*i:2*i+2]
1005      ki2 = 2*ki; //k*6 + i*2
1006      normal[0] = ((double *) normals -> data)[ki2];
1007      normal[1] = ((double *) normals -> data)[ki2+1];
1008      //Edge flux computation
1009      flux_function(ql, qr, zl, zr,
1010                    normal[0], normal[1],
1011                    epsilon, g,
1012                    edgeflux, &max_speed);
1013      //update triangle k
1015      ((double *) stage_explicit_update -> data)[k] -= edgeflux[0]*((double *) edgelengths -> data)[ki];
1016      ((double *) xmom_explicit_update -> data)[k] -= edgeflux[1]*((double *) edgelengths -> data)[ki];
1017      ((double *) ymom_explicit_update -> data)[k] -= edgeflux[2]*((double *) edgelengths -> data)[ki];
1018      //update the neighbour n
1019      if (n>=0){
1021        ((double *) stage_explicit_update -> data)[n] += edgeflux[0]*((double *) edgelengths -> data)[nm];
1022        ((double *) xmom_explicit_update -> data)[n] += edgeflux[1]*((double *) edgelengths -> data)[nm];
1023        ((double *) ymom_explicit_update -> data)[n] += edgeflux[2]*((double *) edgelengths -> data)[nm];
1024      }
1025      ///for (j=0; j<3; j++) {
1026        ///flux[j] -= edgeflux[j]*((double *) edgelengths -> data)[ki];
1027        ///}
1028        //Update timestep
1030        //FIXME: SR Add parameter for CFL condition
1031        if (max_speed > epsilon) {
1032          timestep = min(timestep, ((double *) radii -> data)[k]/max_speed);
1033          //maxspeed in flux_function is calculated as max(|u+a|,|u-a|)
1034          if (n>=0)
1035            timestep = min(timestep, ((double *) radii -> data)[n]/max_speed);
1036        }
1037    } // end for i
1038    //Normalise by area and store for when all conserved
1039    //quantities get updated
1040    ((double *) stage_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
1041    ((double *) xmom_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
1042    ((double *) ymom_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
1043  } //end for k
1044  return Py_BuildValue("d", timestep);
1045}
1046
1047PyObject *protect(PyObject *self, PyObject *args) {
1048  //
1049  //    protect(minimum_allowed_height, maximum_allowed_speed, wc, zc, xmomc, ymomc)
1050
1051
1052  PyArrayObject
1053  *wc,            //Stage at centroids
1054  *zc,            //Elevation at centroids
1055  *xmomc,         //Momentums at centroids
1056  *ymomc;
1057
1058
1059  int N;
1060  double minimum_allowed_height, maximum_allowed_speed, epsilon;
1061
1062  // Convert Python arguments to C
1063  if (!PyArg_ParseTuple(args, "dddOOOO",
1064                        &minimum_allowed_height,
1065                        &maximum_allowed_speed,
1066                        &epsilon,
1067                        &wc, &zc, &xmomc, &ymomc))
1068    return NULL;
1069
1070  N = wc -> dimensions[0];
1071
1072  _protect(N,
1073           minimum_allowed_height,
1074           maximum_allowed_speed,
1075           epsilon,
1076           (double*) wc -> data,
1077           (double*) zc -> data,
1078           (double*) xmomc -> data,
1079           (double*) ymomc -> data);
1080
1081  return Py_BuildValue("");
1082}
1083
1084
1085
1086PyObject *balance_deep_and_shallow(PyObject *self, PyObject *args) {
1087  //
1088  //    balance_deep_and_shallow(wc, zc, hc, wv, zv, hv,
1089  //                             xmomc, ymomc, xmomv, ymomv)
1090
1091
1092  PyArrayObject
1093    *wc,            //Stage at centroids
1094    *zc,            //Elevation at centroids
1095    *hc,            //Height at centroids
1096    *wv,            //Stage at vertices
1097    *zv,            //Elevation at vertices
1098    *hv,            //Depths at vertices
1099    *hvbar,         //h-Limited depths at vertices
1100    *xmomc,         //Momentums at centroids and vertices
1101    *ymomc,
1102    *xmomv,
1103    *ymomv;
1104
1105  int N; //, err;
1106
1107  // Convert Python arguments to C
1108  if (!PyArg_ParseTuple(args, "OOOOOOOOOOO",
1109                        &wc, &zc, &hc,
1110                        &wv, &zv, &hv, &hvbar,
1111                        &xmomc, &ymomc, &xmomv, &ymomv))
1112    return NULL;
1113
1114  N = wc -> dimensions[0];
1115
1116  _balance_deep_and_shallow(N,
1117                            (double*) wc -> data,
1118                            (double*) zc -> data,
1119                            (double*) hc -> data,
1120                            (double*) wv -> data,
1121                            (double*) zv -> data,
1122                            (double*) hv -> data,
1123                            (double*) hvbar -> data,
1124                            (double*) xmomc -> data,
1125                            (double*) ymomc -> data,
1126                            (double*) xmomv -> data,
1127                            (double*) ymomv -> data);
1128
1129
1130  return Py_BuildValue("");
1131}
1132
1133
1134
1135PyObject *h_limiter(PyObject *self, PyObject *args) {
1136
1137  PyObject *domain, *Tmp;
1138  PyArrayObject
1139    *hv, *hc, //Depth at vertices and centroids
1140    *hvbar,   //Limited depth at vertices (return values)
1141    *neighbours;
1142
1143  int k, i, n, N, k3;
1144  int dimensions[2];
1145  double beta_h; //Safety factor (see config.py)
1146  double *hmin, *hmax, hn;
1147
1148  // Convert Python arguments to C
1149  if (!PyArg_ParseTuple(args, "OOO", &domain, &hc, &hv))
1150    return NULL;
1151
1152  neighbours = get_consecutive_array(domain, "neighbours");
1153
1154  //Get safety factor beta_h
1155  Tmp = PyObject_GetAttrString(domain, "beta_h");
1156  if (!Tmp)
1157    return NULL;
1158
1159  beta_h = PyFloat_AsDouble(Tmp);
1160
1161  Py_DECREF(Tmp);
1162
1163  N = hc -> dimensions[0];
1164
1165  //Create hvbar
1166  dimensions[0] = N;
1167  dimensions[1] = 3;
1168  hvbar = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
1169
1170
1171  //Find min and max of this and neighbour's centroid values
1172  hmin = malloc(N * sizeof(double));
1173  hmax = malloc(N * sizeof(double));
1174  for (k=0; k<N; k++) {
1175    k3=k*3;
1176
1177    hmin[k] = ((double*) hc -> data)[k];
1178    hmax[k] = hmin[k];
1179
1180    for (i=0; i<3; i++) {
1181      n = ((long*) neighbours -> data)[k3+i];
1182
1183      //Initialise hvbar with values from hv
1184      ((double*) hvbar -> data)[k3+i] = ((double*) hv -> data)[k3+i];
1185
1186      if (n >= 0) {
1187        hn = ((double*) hc -> data)[n]; //Neighbour's centroid value
1188
1189        hmin[k] = min(hmin[k], hn);
1190        hmax[k] = max(hmax[k], hn);
1191      }
1192    }
1193  }
1194
1195  // Call underlying standard routine
1196  _limit(N, beta_h, (double*) hc -> data, (double*) hvbar -> data, hmin, hmax);
1197
1198  // // //Py_DECREF(domain); //FIXME: NEcessary?
1199  free(hmin);
1200  free(hmax);
1201
1202  //return result using PyArray to avoid memory leak
1203  return PyArray_Return(hvbar);
1204  //return Py_BuildValue("");
1205}
1206
1207PyObject *h_limiter_sw(PyObject *self, PyObject *args) {
1208  //a faster version of h_limiter above
1209  PyObject *domain, *Tmp;
1210  PyArrayObject
1211    *hv, *hc, //Depth at vertices and centroids
1212    *hvbar,   //Limited depth at vertices (return values)
1213    *neighbours;
1214
1215  int k, i, N, k3,k0,k1,k2;
1216  int dimensions[2];
1217  double beta_h; //Safety factor (see config.py)
1218  double hmin, hmax, dh[3];
1219  // Convert Python arguments to C
1220  if (!PyArg_ParseTuple(args, "OOO", &domain, &hc, &hv))
1221    return NULL;
1222  neighbours = get_consecutive_array(domain, "neighbours");
1223
1224  //Get safety factor beta_h
1225  Tmp = PyObject_GetAttrString(domain, "beta_h");
1226  if (!Tmp)
1227    return NULL;
1228  beta_h = PyFloat_AsDouble(Tmp);
1229
1230  Py_DECREF(Tmp);
1231
1232  N = hc -> dimensions[0];
1233
1234  //Create hvbar
1235  dimensions[0] = N;
1236  dimensions[1] = 3;
1237  hvbar = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
1238  for (k=0;k<N;k++){
1239    k3=k*3;
1240    //get the ids of the neighbours
1241    k0 = ((long*) neighbours -> data)[k3];
1242    k1 = ((long*) neighbours -> data)[k3+1];
1243    k2 = ((long*) neighbours -> data)[k3+2];
1244    //set hvbar provisionally
1245    for (i=0;i<3;i++){
1246      ((double*) hvbar -> data)[k3+i] = ((double*) hv -> data)[k3+i];
1247      dh[i]=((double*) hvbar -> data)[k3+i]-((double*) hc -> data)[k];
1248    }
1249    hmin=((double*) hc -> data)[k];
1250    hmax=hmin;
1251    if (k0>=0){
1252      hmin=min(hmin,((double*) hc -> data)[k0]);
1253      hmax=max(hmax,((double*) hc -> data)[k0]);
1254    }
1255    if (k1>=0){
1256      hmin=min(hmin,((double*) hc -> data)[k1]);
1257      hmax=max(hmax,((double*) hc -> data)[k1]);
1258    }
1259    if (k2>=0){
1260      hmin=min(hmin,((double*) hc -> data)[k2]);
1261      hmax=max(hmax,((double*) hc -> data)[k2]);
1262    }
1263    hmin-=((double*) hc -> data)[k];
1264    hmax-=((double*) hc -> data)[k];
1266    for (i=0;i<3;i++)
1267      ((double*) hvbar -> data)[k3+i] = ((double*) hc -> data)[k]+dh[i];
1268  }
1269  return PyArray_Return(hvbar);
1270}
1271
1272PyObject *assign_windfield_values(PyObject *self, PyObject *args) {
1273  //
1274  //      assign_windfield_values(xmom_update, ymom_update,
1275  //                              s_vec, phi_vec, self.const)
1276
1277
1278
1279  PyArrayObject   //(one element per triangle)
1280  *s_vec,         //Speeds
1281  *phi_vec,       //Bearings
1283  *ymom_update;
1284
1285
1286  int N;
1287  double cw;
1288
1289  // Convert Python arguments to C
1290  if (!PyArg_ParseTuple(args, "OOOOd",
1291                        &xmom_update,
1292                        &ymom_update,
1293                        &s_vec, &phi_vec,
1294                        &cw))
1295    return NULL;
1296
1297  N = xmom_update -> dimensions[0];
1298
1299  _assign_wind_field_values(N,
1300           (double*) xmom_update -> data,
1301           (double*) ymom_update -> data,
1302           (double*) s_vec -> data,
1303           (double*) phi_vec -> data,
1304           cw);
1305
1306  return Py_BuildValue("");
1307}
1308
1309
1310
1311
1312//////////////////////////////////////////
1313// Method table for python module
1314static struct PyMethodDef MethodTable[] = {
1315  /* The cast of the function is necessary since PyCFunction values
1316   * only take two PyObject* parameters, and rotate() takes
1317   * three.
1318   */
1319
1320  {"rotate", (PyCFunction)rotate, METH_VARARGS | METH_KEYWORDS, "Print out"},
1321  {"extrapolate_second_order_sw", extrapolate_second_order_sw, METH_VARARGS, "Print out"},
1322  {"compute_fluxes", compute_fluxes, METH_VARARGS, "Print out"},
1323  {"gravity", gravity, METH_VARARGS, "Print out"},
1324  {"manning_friction", manning_friction, METH_VARARGS, "Print out"},
1325  {"balance_deep_and_shallow", balance_deep_and_shallow,
1326   METH_VARARGS, "Print out"},
1327  {"h_limiter", h_limiter,
1328   METH_VARARGS, "Print out"},
1329  {"h_limiter_sw", h_limiter_sw,
1330   METH_VARARGS, "Print out"},
1331  {"protect", protect, METH_VARARGS | METH_KEYWORDS, "Print out"},
1332  {"assign_windfield_values", assign_windfield_values,
1333   METH_VARARGS | METH_KEYWORDS, "Print out"},
1334  //{"distribute_to_vertices_and_edges",
1335  // distribute_to_vertices_and_edges, METH_VARARGS},
1336  //{"update_conserved_quantities",
1337  // update_conserved_quantities, METH_VARARGS},
1338  //{"set_initialcondition",
1339  // set_initialcondition, METH_VARARGS},
1340  {NULL, NULL}
1341};
1342
1343// Module initialisation
1344void initshallow_water_ext(void){
1345  Py_InitModule("shallow_water_ext", MethodTable);
1346
1347  import_array();     //Necessary for handling of NumPY structures
1348}
Note: See TracBrowser for help on using the repository browser.