source: anuga_core/source/anuga/shallow_water/shallow_water_kinetic_ext.c @ 3698

Last change on this file since 3698 was 3698, checked in by steve, 17 years ago

Moved shallow_water_kinetic to shallow_water directory

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