source: inundation/pyvolution/shallow_water_ext.c @ 2380

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