source: anuga_core/source/anuga/shallow_water/shallow_water_ext.c @ 3876

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

Added parameter alpha_balance to config.py to deal with large jumps in bed

File size: 60.6 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
357
358
359void _manning_friction(double g, double eps, int N,
360                       double* w, double* z,
361                       double* uh, double* vh,
362                       double* eta, double* xmom, double* ymom) {
363
364  int k;
365  double S, h;
366
367  for (k=0; k<N; k++) {
368    if (eta[k] > eps) {
369      h = w[k]-z[k];
370      if (h >= eps) {
371        S = -g * eta[k]*eta[k] * sqrt((uh[k]*uh[k] + vh[k]*vh[k]));
372        S /= pow(h, 7.0/3);      //Expensive (on Ole's home computer)
373        //S /= exp(7.0/3.0*log(h));      //seems to save about 15% over manning_friction
374        //S /= h*h*(1 + h/3.0 - h*h/9.0); //FIXME: Could use a Taylor expansion
375
376
377        //Update momentum
378        xmom[k] += S*uh[k];
379        ymom[k] += S*vh[k];
380      }
381    }
382  }
383}
384
385
386/*
387void _manning_friction_explicit(double g, double eps, int N,
388                       double* w, double* z,
389                       double* uh, double* vh,
390                       double* eta, double* xmom, double* ymom) {
391
392  int k;
393  double S, h;
394
395  for (k=0; k<N; k++) {
396    if (eta[k] > eps) {
397      h = w[k]-z[k];
398      if (h >= eps) {
399        S = -g * eta[k]*eta[k] * sqrt((uh[k]*uh[k] + vh[k]*vh[k]));
400        S /= pow(h, 7.0/3);      //Expensive (on Ole's home computer)
401        //S /= exp(7.0/3.0*log(h));      //seems to save about 15% over manning_friction
402        //S /= h*h*(1 + h/3.0 - h*h/9.0); //FIXME: Could use a Taylor expansion
403
404
405        //Update momentum
406        xmom[k] += S*uh[k];
407        ymom[k] += S*vh[k];
408      }
409    }
410  }
411}
412*/
413
414int _balance_deep_and_shallow(int N,
415                              double* wc,
416                              double* zc,
417                              double* hc,
418                              double* wv,
419                              double* zv,
420                              double* hv,
421                              double* hvbar,
422                              double* xmomc,
423                              double* ymomc,
424                              double* xmomv,
425                              double* ymomv,
426                              double alpha_balance) {
427
428  int k, k3, i;
429  double dz, hmin, alpha;
430
431  //Compute linear combination between w-limited stages and
432  //h-limited stages close to the bed elevation.
433
434  for (k=0; k<N; k++) {
435    // Compute maximal variation in bed elevation
436    // This quantitiy is
437    //     dz = max_i abs(z_i - z_c)
438    // and it is independent of dimension
439    // In the 1d case zc = (z0+z1)/2
440    // In the 2d case zc = (z0+z1+z2)/3
441
442    k3 = 3*k;
443
444    //FIXME: Try with this one precomputed
445    dz = 0.0;
446    hmin = hv[k3];
447    for (i=0; i<3; i++) {
448      dz = max(dz, fabs(zv[k3+i]-zc[k]));
449      hmin = min(hmin, hv[k3+i]);
450    }
451
452
453    //Create alpha in [0,1], where alpha==0 means using the h-limited
454    //stage and alpha==1 means using the w-limited stage as
455    //computed by the gradient limiter (both 1st or 2nd order)
456    //
457    //If hmin > dz/2 then alpha = 1 and the bed will have no effect
458    //If hmin < 0 then alpha = 0 reverting to constant height above bed.
459
460
461    if (dz > 0.0)
462      //if (hmin<0.0)
463      //        alpha = 0.0;
464      //else
465      //  alpha = max( min( hc[k]/dz, 1.0), 0.0 );
466      alpha = max( min( alpha_balance*hmin/dz, 1.0), 0.0 );
467    else
468      alpha = 1.0;  //Flat bed
469     
470     
471    //alpha = 1.0;  //Always deep FIXME: This actually looks good now
472
473    //printf("dz = %.3f, alpha = %.8f\n", dz, alpha);
474
475    //  Let
476    //
477    //    wvi be the w-limited stage (wvi = zvi + hvi)
478    //    wvi- be the h-limited state (wvi- = zvi + hvi-)
479    //
480    //
481    //  where i=0,1,2 denotes the vertex ids
482    //
483    //  Weighted balance between w-limited and h-limited stage is
484    //
485    //    wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi)
486    //
487    //  It follows that the updated wvi is
488    //    wvi := zvi + (1-alpha)*hvi- + alpha*hvi
489    //
490    //   Momentum is balanced between constant and limited
491
492    if (alpha < 1) {
493      for (i=0; i<3; i++) {
494         wv[k3+i] = zv[k3+i] + (1-alpha)*hvbar[k3+i] + alpha*hv[k3+i];
495
496        //Update momentum as a linear combination of
497        //xmomc and ymomc (shallow) and momentum
498        //from extrapolator xmomv and ymomv (deep).
499        xmomv[k3+i] = (1-alpha)*xmomc[k] + alpha*xmomv[k3+i];
500        ymomv[k3+i] = (1-alpha)*ymomc[k] + alpha*ymomv[k3+i];
501      }
502    }
503  }
504  return 0;
505}
506
507
508
509int _protect(int N,
510             double minimum_allowed_height,
511             double maximum_allowed_speed,
512             double epsilon,
513             double* wc,
514             double* zc,
515             double* xmomc,
516             double* ymomc) {
517
518  int k;
519  double hc;
520  double u, v, reduced_speed;
521
522  //Protect against initesimal and negative heights
523  for (k=0; k<N; k++) {
524    hc = wc[k] - zc[k];
525
526    if (hc < minimum_allowed_height) {
527       
528      //Old code: Set momentum to zero and ensure h is non negative
529      //xmomc[k] = 0.0;
530      //ymomc[k] = 0.0;
531      //if (hc <= 0.0) wc[k] = zc[k];
532
533
534      //New code: Adjust momentum to guarantee speeds are physical
535      //          ensure h is non negative
536      //FIXME (Ole): This is only implemented in this C extension and
537      //             has no Python equivalent
538           
539      if (hc <= 0.0) {
540        wc[k] = zc[k];
541        xmomc[k] = 0.0;
542        ymomc[k] = 0.0;
543      } else {
544        //Reduce excessive speeds derived from division by small hc
545       
546        u = xmomc[k]/hc;
547        if (fabs(u) > maximum_allowed_speed) {
548          reduced_speed = maximum_allowed_speed * u/fabs(u);
549          //printf("Speed (u) has been reduced from %.3f to %.3f\n",
550          //     u, reduced_speed);
551          xmomc[k] = reduced_speed * hc;
552        }
553
554        v = ymomc[k]/hc;
555        if (fabs(v) > maximum_allowed_speed) {
556          reduced_speed = maximum_allowed_speed * v/fabs(v);
557          //printf("Speed (v) has been reduced from %.3f to %.3f\n",
558          //     v, reduced_speed);
559          ymomc[k] = reduced_speed * hc;
560        }
561      }
562    }
563  }
564  return 0;
565}
566
567
568
569
570int _assign_wind_field_values(int N,
571                              double* xmom_update,
572                              double* ymom_update,
573                              double* s_vec,
574                              double* phi_vec,
575                              double cw) {
576
577  //Assign windfield values to momentum updates
578
579  int k;
580  double S, s, phi, u, v;
581
582  for (k=0; k<N; k++) {
583
584    s = s_vec[k];
585    phi = phi_vec[k];
586
587    //Convert to radians
588    phi = phi*pi/180;
589
590    //Compute velocity vector (u, v)
591    u = s*cos(phi);
592    v = s*sin(phi);
593
594    //Compute wind stress
595    S = cw * sqrt(u*u + v*v);
596    xmom_update[k] += S*u;
597    ymom_update[k] += S*v;
598  }
599  return 0;
600}
601
602
603
604///////////////////////////////////////////////////////////////////
605// Gateways to Python
606
607PyObject *gravity(PyObject *self, PyObject *args) {
608  //
609  //  gravity(g, h, v, x, xmom, ymom)
610  //
611
612
613  PyArrayObject *h, *v, *x, *xmom, *ymom;
614  int k, i, N, k3, k6;
615  double g, avg_h, zx, zy;
616  double x0, y0, x1, y1, x2, y2, z0, z1, z2;
617
618  if (!PyArg_ParseTuple(args, "dOOOOO",
619                        &g, &h, &v, &x,
620                        &xmom, &ymom)) {
621    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: gravity could not parse input arguments");
622    return NULL;
623  }
624
625  N = h -> dimensions[0];
626  for (k=0; k<N; k++) {
627    k3 = 3*k;  // base index
628    k6 = 6*k;  // base index
629
630    avg_h = 0.0;
631    for (i=0; i<3; i++) {
632      avg_h += ((double *) h -> data)[k3+i];
633    }
634    avg_h /= 3;
635
636
637    //Compute bed slope
638    x0 = ((double*) x -> data)[k6 + 0];
639    y0 = ((double*) x -> data)[k6 + 1];
640    x1 = ((double*) x -> data)[k6 + 2];
641    y1 = ((double*) x -> data)[k6 + 3];
642    x2 = ((double*) x -> data)[k6 + 4];
643    y2 = ((double*) x -> data)[k6 + 5];
644
645
646    z0 = ((double*) v -> data)[k3 + 0];
647    z1 = ((double*) v -> data)[k3 + 1];
648    z2 = ((double*) v -> data)[k3 + 2];
649
650    _gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2, &zx, &zy);
651
652    //Update momentum
653    ((double*) xmom -> data)[k] += -g*zx*avg_h;
654    ((double*) ymom -> data)[k] += -g*zy*avg_h;
655  }
656
657  return Py_BuildValue("");
658}
659
660
661PyObject *manning_friction(PyObject *self, PyObject *args) {
662  //
663  // manning_friction(g, eps, h, uh, vh, eta, xmom_update, ymom_update)
664  //
665
666
667  PyArrayObject *w, *z, *uh, *vh, *eta, *xmom, *ymom;
668  int N;
669  double g, eps;
670
671  if (!PyArg_ParseTuple(args, "ddOOOOOOO",
672                        &g, &eps, &w, &z, &uh, &vh, &eta,
673                        &xmom, &ymom)) {
674    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: manning_friction could not parse input arguments");
675    return NULL;
676  }
677
678
679  N = w -> dimensions[0];
680  _manning_friction(g, eps, N,
681                    (double*) w -> data,
682                    (double*) z -> data,
683                    (double*) uh -> data,
684                    (double*) vh -> data,
685                    (double*) eta -> data,
686                    (double*) xmom -> data,
687                    (double*) ymom -> data);
688
689  return Py_BuildValue("");
690}
691
692
693/*
694PyObject *manning_friction_explicit(PyObject *self, PyObject *args) {
695  //
696  // manning_friction_explicit(g, eps, h, uh, vh, eta, xmom_update, ymom_update)
697  //
698
699
700  PyArrayObject *w, *z, *uh, *vh, *eta, *xmom, *ymom;
701  int N;
702  double g, eps;
703
704  if (!PyArg_ParseTuple(args, "ddOOOOOOO",
705                        &g, &eps, &w, &z, &uh, &vh, &eta,
706                        &xmom, &ymom))
707    return NULL;
708
709  N = w -> dimensions[0];
710  _manning_friction_explicit(g, eps, N,
711                    (double*) w -> data,
712                    (double*) z -> data,
713                    (double*) uh -> data,
714                    (double*) vh -> data,
715                    (double*) eta -> data,
716                    (double*) xmom -> data,
717                    (double*) ymom -> data);
718
719  return Py_BuildValue("");
720}
721*/
722
723PyObject *extrapolate_second_order_sw(PyObject *self, PyObject *args) {
724  /*Compute the vertex values based on a linear reconstruction on each triangle
725    These values are calculated as follows:
726    1) For each triangle not adjacent to a boundary, we consider the auxiliary triangle
727    formed by the centroids of its three neighbours.
728    2) For each conserved quantity, we integrate around the auxiliary triangle's boundary the product
729    of the quantity and the outward normal vector. Dividing by the triangle area gives (a,b), the average
730    of the vector (q_x,q_y) on the auxiliary triangle. We suppose that the linear reconstruction on the
731    original triangle has gradient (a,b).
732    3) Provisional vertex junmps dqv[0,1,2] are computed and these are then limited by calling the functions
733    find_qmin_and_qmax and limit_gradient
734
735    Python call:
736    extrapolate_second_order_sw(domain.surrogate_neighbours,
737                                domain.number_of_boundaries
738                                domain.centroid_coordinates,
739                                Stage.centroid_values
740                                Xmom.centroid_values
741                                Ymom.centroid_values
742                                domain.vertex_coordinates,
743                                Stage.vertex_values,
744                                Xmom.vertex_values,
745                                Ymom.vertex_values)
746
747    Post conditions:
748            The vertices of each triangle have values from a limited linear reconstruction
749            based on centroid values
750
751  */
752  PyArrayObject *surrogate_neighbours,
753    *number_of_boundaries,
754    *centroid_coordinates,
755    *stage_centroid_values,
756    *xmom_centroid_values,
757    *ymom_centroid_values,
758        *elevation_centroid_values,
759    *vertex_coordinates,
760    *stage_vertex_values,
761    *xmom_vertex_values,
762    *ymom_vertex_values,
763        *elevation_vertex_values;
764  PyObject *domain, *Tmp;
765  //Local variables
766  double a, b;//gradient vector, not stored but used to calculate vertex values from centroids
767  int number_of_elements,k,k0,k1,k2,k3,k6,coord_index,i;
768  double x,y,x0,y0,x1,y1,x2,y2,xv0,yv0,xv1,yv1,xv2,yv2;//vertices of the auxiliary triangle
769  double dx1,dx2,dy1,dy2,dxv0,dxv1,dxv2,dyv0,dyv1,dyv2,dq0,dq1,dq2,area2;
770  double dqv[3], qmin, qmax, hmin;
771  double hc, h0, h1, h2;
772  double beta_w, beta_w_dry, beta_uh, beta_uh_dry, beta_vh, beta_vh_dry, beta_tmp;
773  double minimum_allowed_height;
774  //provisional jumps from centroids to v'tices and safety factor re limiting
775  //by which these jumps are limited
776  // Convert Python arguments to C
777  if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOO",
778                        &domain,
779                        &surrogate_neighbours,
780                        &number_of_boundaries,
781                        &centroid_coordinates,
782                        &stage_centroid_values,
783                        &xmom_centroid_values,
784                        &ymom_centroid_values,
785                        &elevation_centroid_values,
786                        &vertex_coordinates,
787                        &stage_vertex_values,
788                        &xmom_vertex_values,
789                        &ymom_vertex_values,
790                        &elevation_vertex_values)) {
791    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
792    return NULL;
793  }
794
795  //get the safety factor beta_w, set in the config.py file. This is used in the limiting process
796  Tmp = PyObject_GetAttrString(domain, "beta_w");
797  if (!Tmp) {
798    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_w from domain");
799    return NULL;
800  } 
801  beta_w = PyFloat_AsDouble(Tmp);
802  Py_DECREF(Tmp);
803 
804  Tmp = PyObject_GetAttrString(domain, "beta_w_dry");
805  if (!Tmp) {
806    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_w_dry from domain");
807    return NULL;
808  } 
809  beta_w_dry = PyFloat_AsDouble(Tmp);
810  Py_DECREF(Tmp);
811 
812  Tmp = PyObject_GetAttrString(domain, "beta_uh");
813  if (!Tmp) {
814    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_uh from domain");
815    return NULL;
816  } 
817  beta_uh = PyFloat_AsDouble(Tmp);
818  Py_DECREF(Tmp);
819 
820  Tmp = PyObject_GetAttrString(domain, "beta_uh_dry");
821  if (!Tmp) {
822    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_uh_dry from domain");
823    return NULL;
824  } 
825  beta_uh_dry = PyFloat_AsDouble(Tmp);
826  Py_DECREF(Tmp); 
827
828  Tmp = PyObject_GetAttrString(domain, "beta_vh");
829  if (!Tmp) {
830    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_vh from domain");
831    return NULL;
832  } 
833  beta_vh = PyFloat_AsDouble(Tmp);
834  Py_DECREF(Tmp);
835 
836  Tmp = PyObject_GetAttrString(domain, "beta_vh_dry");
837  if (!Tmp) {
838    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_vh_dry from domain");
839    return NULL;
840  } 
841  beta_vh_dry = PyFloat_AsDouble(Tmp);
842  Py_DECREF(Tmp);
843 
844  Tmp = PyObject_GetAttrString(domain, "minimum_allowed_height");
845  if (!Tmp) {
846    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object minimum_allowed_heigt");
847    return NULL;
848  } 
849  minimum_allowed_height = PyFloat_AsDouble(Tmp);
850  Py_DECREF(Tmp); 
851 
852  number_of_elements = stage_centroid_values -> dimensions[0];
853  for (k=0; k<number_of_elements; k++) {
854    k3=k*3;
855    k6=k*6;
856
857    if (((long *) number_of_boundaries->data)[k]==3){/*no neighbours, set gradient on the triangle to zero*/
858      ((double *) stage_vertex_values->data)[k3]=((double *)stage_centroid_values->data)[k];
859      ((double *) stage_vertex_values->data)[k3+1]=((double *)stage_centroid_values->data)[k];
860      ((double *) stage_vertex_values->data)[k3+2]=((double *)stage_centroid_values->data)[k];
861      ((double *) xmom_vertex_values->data)[k3]=((double *)xmom_centroid_values->data)[k];
862      ((double *) xmom_vertex_values->data)[k3+1]=((double *)xmom_centroid_values->data)[k];
863      ((double *) xmom_vertex_values->data)[k3+2]=((double *)xmom_centroid_values->data)[k];
864      ((double *) ymom_vertex_values->data)[k3]=((double *)ymom_centroid_values->data)[k];
865      ((double *) ymom_vertex_values->data)[k3+1]=((double *)ymom_centroid_values->data)[k];
866      ((double *) ymom_vertex_values->data)[k3+2]=((double *)ymom_centroid_values->data)[k];
867      continue;
868    }
869    else {//we will need centroid coordinates and vertex coordinates of the triangle
870      //get the vertex coordinates of the FV triangle
871      xv0=((double *)vertex_coordinates->data)[k6]; yv0=((double *)vertex_coordinates->data)[k6+1];
872      xv1=((double *)vertex_coordinates->data)[k6+2]; yv1=((double *)vertex_coordinates->data)[k6+3];
873      xv2=((double *)vertex_coordinates->data)[k6+4]; yv2=((double *)vertex_coordinates->data)[k6+5];
874      //get the centroid coordinates of the FV triangle
875      coord_index=2*k;
876      x=((double *)centroid_coordinates->data)[coord_index];
877      y=((double *)centroid_coordinates->data)[coord_index+1];
878      //store x- and y- differentials for the vertices of the FV triangle relative to the centroid
879      dxv0=xv0-x; dxv1=xv1-x; dxv2=xv2-x;
880      dyv0=yv0-y; dyv1=yv1-y; dyv2=yv2-y;
881    }
882    if (((long *)number_of_boundaries->data)[k]<=1){
883      //if no boundaries, auxiliary triangle is formed from the centroids of the three neighbours
884      //if one boundary, auxiliary triangle is formed from this centroid and its two neighbours
885      k0=((long *)surrogate_neighbours->data)[k3];
886      k1=((long *)surrogate_neighbours->data)[k3+1];
887      k2=((long *)surrogate_neighbours->data)[k3+2];
888      //get the auxiliary triangle's vertex coordinates (really the centroids of neighbouring triangles)
889      coord_index=2*k0;
890      x0=((double *)centroid_coordinates->data)[coord_index];
891      y0=((double *)centroid_coordinates->data)[coord_index+1];
892      coord_index=2*k1;
893      x1=((double *)centroid_coordinates->data)[coord_index];
894      y1=((double *)centroid_coordinates->data)[coord_index+1];
895      coord_index=2*k2;
896      x2=((double *)centroid_coordinates->data)[coord_index];
897      y2=((double *)centroid_coordinates->data)[coord_index+1];
898      //store x- and y- differentials for the vertices of the auxiliary triangle
899      dx1=x1-x0; dx2=x2-x0;
900      dy1=y1-y0; dy2=y2-y0;
901      //calculate 2*area of the auxiliary triangle
902      area2 = dy2*dx1 - dy1*dx2;//the triangle is guaranteed to be counter-clockwise
903      //If the mesh is 'weird' near the boundary, the trianlge might be flat or clockwise:
904      if (area2<=0) {
905        PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: negative triangle area encountered");
906        return NULL;
907      } 
908     
909
910      //### Calculate heights of neighbouring cells
911      hc = ((double *)stage_centroid_values->data)[k]  - ((double *)elevation_centroid_values->data)[k];
912      h0 = ((double *)stage_centroid_values->data)[k0] - ((double *)elevation_centroid_values->data)[k0];
913      h1 = ((double *)stage_centroid_values->data)[k1] - ((double *)elevation_centroid_values->data)[k1];
914      h2 = ((double *)stage_centroid_values->data)[k2] - ((double *)elevation_centroid_values->data)[k2];
915      hmin = min(hc,min(h0,min(h1,h2)));
916     
917      //### stage ###
918      //calculate the difference between vertex 0 of the auxiliary triangle and the FV triangle centroid
919      dq0=((double *)stage_centroid_values->data)[k0]-((double *)stage_centroid_values->data)[k];
920      //calculate differentials between the vertices of the auxiliary triangle
921      dq1=((double *)stage_centroid_values->data)[k1]-((double *)stage_centroid_values->data)[k0];
922      dq2=((double *)stage_centroid_values->data)[k2]-((double *)stage_centroid_values->data)[k0];
923      //calculate the gradient of stage on the auxiliary triangle
924      a = dy2*dq1 - dy1*dq2;
925      a /= area2;
926      b = dx1*dq2 - dx2*dq1;
927      b /= area2;
928      //calculate provisional jumps in stage from the centroid of the FV tri to its vertices, to be limited
929      dqv[0]=a*dxv0+b*dyv0;
930      dqv[1]=a*dxv1+b*dyv1;
931      dqv[2]=a*dxv2+b*dyv2;
932      //now we want to find min and max of the centroid and the vertices of the auxiliary triangle
933      //and compute jumps from the centroid to the min and max
934      find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
935      // Playing with dry wet interface
936      hmin = qmin;
937      beta_tmp = beta_w;
938      if (hmin<minimum_allowed_height)
939        beta_tmp = beta_w_dry;
940      limit_gradient(dqv,qmin,qmax,beta_tmp);//the gradient will be limited
941      for (i=0;i<3;i++)
942        ((double *)stage_vertex_values->data)[k3+i]=((double *)stage_centroid_values->data)[k]+dqv[i];
943     
944      //### xmom ###
945      //calculate the difference between vertex 0 of the auxiliary triangle and the FV triangle centroid
946      dq0=((double *)xmom_centroid_values->data)[k0]-((double *)xmom_centroid_values->data)[k];
947      //calculate differentials between the vertices of the auxiliary triangle
948      dq1=((double *)xmom_centroid_values->data)[k1]-((double *)xmom_centroid_values->data)[k0];
949      dq2=((double *)xmom_centroid_values->data)[k2]-((double *)xmom_centroid_values->data)[k0];
950      //calculate the gradient of xmom on the auxiliary triangle
951      a = dy2*dq1 - dy1*dq2;
952      a /= area2;
953      b = dx1*dq2 - dx2*dq1;
954      b /= area2;
955      //calculate provisional jumps in stage from the centroid of the FV tri to its vertices, to be limited
956      dqv[0]=a*dxv0+b*dyv0;
957      dqv[1]=a*dxv1+b*dyv1;
958      dqv[2]=a*dxv2+b*dyv2;
959      //now we want to find min and max of the centroid and the vertices of the auxiliary triangle
960      //and compute jumps from the centroid to the min and max
961      find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
962      beta_tmp = beta_uh;
963      if (hmin<minimum_allowed_height)
964        beta_tmp = beta_uh_dry;
965      limit_gradient(dqv,qmin,qmax,beta_tmp);//the gradient will be limited
966      for (i=0;i<3;i++)
967        ((double *)xmom_vertex_values->data)[k3+i]=((double *)xmom_centroid_values->data)[k]+dqv[i];
968     
969      //### ymom ###
970      //calculate the difference between vertex 0 of the auxiliary triangle and the FV triangle centroid
971      dq0=((double *)ymom_centroid_values->data)[k0]-((double *)ymom_centroid_values->data)[k];
972      //calculate differentials between the vertices of the auxiliary triangle
973      dq1=((double *)ymom_centroid_values->data)[k1]-((double *)ymom_centroid_values->data)[k0];
974      dq2=((double *)ymom_centroid_values->data)[k2]-((double *)ymom_centroid_values->data)[k0];
975      //calculate the gradient of xmom on the auxiliary triangle
976      a = dy2*dq1 - dy1*dq2;
977      a /= area2;
978      b = dx1*dq2 - dx2*dq1;
979      b /= area2;
980      //calculate provisional jumps in stage from the centroid of the FV tri to its vertices, to be limited
981      dqv[0]=a*dxv0+b*dyv0;
982      dqv[1]=a*dxv1+b*dyv1;
983      dqv[2]=a*dxv2+b*dyv2;
984      //now we want to find min and max of the centroid and the vertices of the auxiliary triangle
985      //and compute jumps from the centroid to the min and max
986      find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
987      beta_tmp = beta_vh;
988      if (hmin<minimum_allowed_height)
989        beta_tmp = beta_vh_dry;
990      limit_gradient(dqv,qmin,qmax,beta_tmp);//the gradient will be limited
991      for (i=0;i<3;i++)
992        ((double *)ymom_vertex_values->data)[k3+i]=((double *)ymom_centroid_values->data)[k]+dqv[i];
993    }//if (number_of_boundaries[k]<=1)
994    else{//number_of_boundaries==2
995      //one internal neighbour and gradient is in direction of the neighbour's centroid
996      //find the only internal neighbour
997      for (k2=k3;k2<k3+3;k2++){//k2 just indexes the edges of triangle k
998        if (((long *)surrogate_neighbours->data)[k2]!=k)//find internal neighbour of triabngle k
999          break;
1000      }
1001      if ((k2==k3+3)) {//if we didn't find an internal neighbour
1002        PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: Internal neighbour not found");     
1003        return NULL;//error
1004      }
1005     
1006      k1=((long *)surrogate_neighbours->data)[k2];
1007      //the coordinates of the triangle are already (x,y). Get centroid of the neighbour (x1,y1)
1008      coord_index=2*k1;
1009      x1=((double *)centroid_coordinates->data)[coord_index];
1010      y1=((double *)centroid_coordinates->data)[coord_index+1];
1011      //compute x- and y- distances between the centroid of the FV triangle and that of its neighbour
1012      dx1=x1-x; dy1=y1-y;
1013      //set area2 as the square of the distance
1014      area2=dx1*dx1+dy1*dy1;
1015      //set dx2=(x1-x0)/((x1-x0)^2+(y1-y0)^2) and dy2=(y1-y0)/((x1-x0)^2+(y1-y0)^2) which
1016      //respectively correspond to the x- and y- gradients of the conserved quantities
1017      dx2=1.0/area2;
1018      dy2=dx2*dy1;
1019      dx2*=dx1;
1020     
1021      //## stage ###
1022      //compute differentials
1023      dq1=((double *)stage_centroid_values->data)[k1]-((double *)stage_centroid_values->data)[k];
1024      //calculate the gradient between the centroid of the FV triangle and that of its neighbour
1025      a=dq1*dx2;
1026      b=dq1*dy2;
1027      //calculate provisional vertex jumps, to be limited
1028      dqv[0]=a*dxv0+b*dyv0;
1029      dqv[1]=a*dxv1+b*dyv1;
1030      dqv[2]=a*dxv2+b*dyv2;
1031      //now limit the jumps
1032      if (dq1>=0.0){
1033        qmin=0.0;
1034        qmax=dq1;
1035      }
1036      else{
1037        qmin=dq1;
1038        qmax=0.0;
1039      }
1040     
1041     
1042      limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
1043      for (i=0;i<3;i++)
1044        ((double *)stage_vertex_values->data)[k3+i]=((double *)stage_centroid_values->data)[k]+dqv[i];
1045     
1046      //## xmom ###
1047      //compute differentials
1048      dq1=((double *)xmom_centroid_values->data)[k1]-((double *)xmom_centroid_values->data)[k];
1049      //calculate the gradient between the centroid of the FV triangle and that of its neighbour
1050      a=dq1*dx2;
1051      b=dq1*dy2;
1052      //calculate provisional vertex jumps, to be limited
1053      dqv[0]=a*dxv0+b*dyv0;
1054      dqv[1]=a*dxv1+b*dyv1;
1055      dqv[2]=a*dxv2+b*dyv2;
1056      //now limit the jumps
1057      if (dq1>=0.0){
1058        qmin=0.0;
1059        qmax=dq1;
1060      }
1061      else{
1062        qmin=dq1;
1063        qmax=0.0;
1064      }
1065      limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
1066      for (i=0;i<3;i++)
1067        ((double *)xmom_vertex_values->data)[k3+i]=((double *)xmom_centroid_values->data)[k]+dqv[i];
1068     
1069      //## ymom ###
1070      //compute differentials
1071      dq1=((double *)ymom_centroid_values->data)[k1]-((double *)ymom_centroid_values->data)[k];
1072      //calculate the gradient between the centroid of the FV triangle and that of its neighbour
1073      a=dq1*dx2;
1074      b=dq1*dy2;
1075      //calculate provisional vertex jumps, to be limited
1076      dqv[0]=a*dxv0+b*dyv0;
1077      dqv[1]=a*dxv1+b*dyv1;
1078      dqv[2]=a*dxv2+b*dyv2;
1079      //now limit the jumps
1080      if (dq1>=0.0){
1081        qmin=0.0;
1082        qmax=dq1;
1083      }
1084      else{
1085        qmin=dq1;
1086        qmax=0.0;
1087      }
1088      limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
1089      for (i=0;i<3;i++)
1090        ((double *)ymom_vertex_values->data)[k3+i]=((double *)ymom_centroid_values->data)[k]+dqv[i];
1091    }//else [number_of_boudaries==2]
1092  }//for k=0 to number_of_elements-1
1093  return Py_BuildValue("");
1094}//extrapolate_second-order_sw
1095
1096
1097PyObject *rotate(PyObject *self, PyObject *args, PyObject *kwargs) {
1098  //
1099  // r = rotate(q, normal, direction=1)
1100  //
1101  // Where q is assumed to be a Float numeric array of length 3 and
1102  // normal a Float numeric array of length 2.
1103
1104
1105  PyObject *Q, *Normal;
1106  PyArrayObject *q, *r, *normal;
1107
1108  static char *argnames[] = {"q", "normal", "direction", NULL};
1109  int dimensions[1], i, direction=1;
1110  double n1, n2;
1111
1112  // Convert Python arguments to C
1113  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OO|i", argnames,
1114                                   &Q, &Normal, &direction)) {
1115    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: rotate could not parse input arguments");
1116    return NULL;
1117  } 
1118
1119  //Input checks (convert sequences into numeric arrays)
1120  q = (PyArrayObject *)
1121    PyArray_ContiguousFromObject(Q, PyArray_DOUBLE, 0, 0);
1122  normal = (PyArrayObject *)
1123    PyArray_ContiguousFromObject(Normal, PyArray_DOUBLE, 0, 0);
1124
1125
1126  if (normal -> dimensions[0] != 2) {
1127    PyErr_SetString(PyExc_RuntimeError, "Normal vector must have 2 components");
1128    return NULL;
1129  }
1130
1131  //Allocate space for return vector r (don't DECREF)
1132  dimensions[0] = 3;
1133  r = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
1134
1135  //Copy
1136  for (i=0; i<3; i++) {
1137    ((double *) (r -> data))[i] = ((double *) (q -> data))[i];
1138  }
1139
1140  //Get normal and direction
1141  n1 = ((double *) normal -> data)[0];
1142  n2 = ((double *) normal -> data)[1];
1143  if (direction == -1) n2 = -n2;
1144
1145  //Rotate
1146  _rotate((double *) r -> data, n1, n2);
1147
1148  //Release numeric arrays
1149  Py_DECREF(q);
1150  Py_DECREF(normal);
1151
1152  //return result using PyArray to avoid memory leak
1153  return PyArray_Return(r);
1154}
1155
1156PyObject *compute_fluxes_ext_central(PyObject *self, PyObject *args) {
1157  /*Compute all fluxes and the timestep suitable for all volumes
1158    in domain.
1159
1160    Compute total flux for each conserved quantity using "flux_function_central"
1161
1162    Fluxes across each edge are scaled by edgelengths and summed up
1163    Resulting flux is then scaled by area and stored in
1164    explicit_update for each of the three conserved quantities
1165    stage, xmomentum and ymomentum
1166
1167    The maximal allowable speed computed by the flux_function for each volume
1168    is converted to a timestep that must not be exceeded. The minimum of
1169    those is computed as the next overall timestep.
1170
1171    Python call:
1172    domain.timestep = compute_fluxes(timestep,
1173                                     domain.epsilon,
1174                                     domain.g,
1175                                     domain.neighbours,
1176                                     domain.neighbour_edges,
1177                                     domain.normals,
1178                                     domain.edgelengths,
1179                                     domain.radii,
1180                                     domain.areas,
1181                                     Stage.edge_values,
1182                                     Xmom.edge_values,
1183                                     Ymom.edge_values,
1184                                     Bed.edge_values,
1185                                     Stage.boundary_values,
1186                                     Xmom.boundary_values,
1187                                     Ymom.boundary_values,
1188                                     Stage.explicit_update,
1189                                     Xmom.explicit_update,
1190                                     Ymom.explicit_update,
1191                                     already_computed_flux)
1192
1193
1194    Post conditions:
1195      domain.explicit_update is reset to computed flux values
1196      domain.timestep is set to the largest step satisfying all volumes.
1197
1198
1199  */
1200
1201
1202  PyArrayObject *neighbours, *neighbour_edges,
1203    *normals, *edgelengths, *radii, *areas,
1204    *tri_full_flag,
1205    *stage_edge_values,
1206    *xmom_edge_values,
1207    *ymom_edge_values,
1208    *bed_edge_values,
1209    *stage_boundary_values,
1210    *xmom_boundary_values,
1211    *ymom_boundary_values,
1212    *stage_explicit_update,
1213    *xmom_explicit_update,
1214    *ymom_explicit_update,
1215    *already_computed_flux;//tracks whether the flux across an edge has already been computed
1216
1217
1218  //Local variables
1219  double timestep, max_speed, epsilon, g;
1220  double normal[2], ql[3], qr[3], zl, zr;
1221  double edgeflux[3]; //Work arrays for summing up fluxes
1222
1223  int number_of_elements, k, i, m, n;
1224  int ki, nm=0, ki2; //Index shorthands
1225  static long call=1;
1226
1227
1228  // Convert Python arguments to C
1229  if (!PyArg_ParseTuple(args, "dddOOOOOOOOOOOOOOOOOO",
1230                        &timestep,
1231                        &epsilon,
1232                        &g,
1233                        &neighbours,
1234                        &neighbour_edges,
1235                        &normals,
1236                        &edgelengths, &radii, &areas,
1237                        &tri_full_flag,
1238                        &stage_edge_values,
1239                        &xmom_edge_values,
1240                        &ymom_edge_values,
1241                        &bed_edge_values,
1242                        &stage_boundary_values,
1243                        &xmom_boundary_values,
1244                        &ymom_boundary_values,
1245                        &stage_explicit_update,
1246                        &xmom_explicit_update,
1247                        &ymom_explicit_update,
1248                        &already_computed_flux)) {
1249    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
1250    return NULL;
1251  }
1252  number_of_elements = stage_edge_values -> dimensions[0];
1253  call++;//a static local variable to which already_computed_flux is compared
1254  //set explicit_update to zero for all conserved_quantities.
1255  //This assumes compute_fluxes called before forcing terms
1256  for (k=0; k<number_of_elements; k++) {
1257    ((double *) stage_explicit_update -> data)[k]=0.0;
1258    ((double *) xmom_explicit_update -> data)[k]=0.0;
1259    ((double *) ymom_explicit_update -> data)[k]=0.0;
1260  }
1261  //Loop through neighbours and compute edge flux for each
1262  for (k=0; k<number_of_elements; k++) {
1263    for (i=0; i<3; i++) {
1264      ki = k*3+i;
1265      if (((long *) already_computed_flux->data)[ki]==call)//we've already computed the flux across this edge
1266        continue;
1267      ql[0] = ((double *) stage_edge_values -> data)[ki];
1268      ql[1] = ((double *) xmom_edge_values -> data)[ki];
1269      ql[2] = ((double *) ymom_edge_values -> data)[ki];
1270      zl =    ((double *) bed_edge_values -> data)[ki];
1271
1272      //Quantities at neighbour on nearest face
1273      n = ((long *) neighbours -> data)[ki];
1274      if (n < 0) {
1275        m = -n-1; //Convert negative flag to index
1276        qr[0] = ((double *) stage_boundary_values -> data)[m];
1277        qr[1] = ((double *) xmom_boundary_values -> data)[m];
1278        qr[2] = ((double *) ymom_boundary_values -> data)[m];
1279        zr = zl; //Extend bed elevation to boundary
1280      } else {
1281        m = ((long *) neighbour_edges -> data)[ki];
1282        nm = n*3+m;
1283        qr[0] = ((double *) stage_edge_values -> data)[nm];
1284        qr[1] = ((double *) xmom_edge_values -> data)[nm];
1285        qr[2] = ((double *) ymom_edge_values -> data)[nm];
1286        zr =    ((double *) bed_edge_values -> data)[nm];
1287      }
1288      // Outward pointing normal vector
1289      // normal = domain.normals[k, 2*i:2*i+2]
1290      ki2 = 2*ki; //k*6 + i*2
1291      normal[0] = ((double *) normals -> data)[ki2];
1292      normal[1] = ((double *) normals -> data)[ki2+1];
1293      //Edge flux computation
1294      flux_function_central(ql, qr, zl, zr,
1295                    normal[0], normal[1],
1296                    epsilon, g,
1297                    edgeflux, &max_speed);
1298      //update triangle k
1299      ((long *) already_computed_flux->data)[ki]=call;
1300      ((double *) stage_explicit_update -> data)[k] -= edgeflux[0]*((double *) edgelengths -> data)[ki];
1301      ((double *) xmom_explicit_update -> data)[k] -= edgeflux[1]*((double *) edgelengths -> data)[ki];
1302      ((double *) ymom_explicit_update -> data)[k] -= edgeflux[2]*((double *) edgelengths -> data)[ki];
1303      //update the neighbour n
1304      if (n>=0){
1305        ((long *) already_computed_flux->data)[nm]=call;
1306        ((double *) stage_explicit_update -> data)[n] += edgeflux[0]*((double *) edgelengths -> data)[nm];
1307        ((double *) xmom_explicit_update -> data)[n] += edgeflux[1]*((double *) edgelengths -> data)[nm];
1308        ((double *) ymom_explicit_update -> data)[n] += edgeflux[2]*((double *) edgelengths -> data)[nm];
1309      }
1310      ///for (j=0; j<3; j++) {
1311        ///flux[j] -= edgeflux[j]*((double *) edgelengths -> data)[ki];
1312        ///}
1313        //Update timestep
1314        //timestep = min(timestep, domain.radii[k]/max_speed)
1315        //FIXME: SR Add parameter for CFL condition
1316    if ( ((long *) tri_full_flag -> data)[k] == 1) {
1317            if (max_speed > epsilon) {
1318                timestep = min(timestep, ((double *) radii -> data)[k]/max_speed);
1319                //maxspeed in flux_function is calculated as max(|u+a|,|u-a|)
1320                if (n>=0)
1321                    timestep = min(timestep, ((double *) radii -> data)[n]/max_speed);
1322            }
1323    }
1324    } // end for i
1325    //Normalise by area and store for when all conserved
1326    //quantities get updated
1327    ((double *) stage_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
1328    ((double *) xmom_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
1329    ((double *) ymom_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
1330  } //end for k
1331  return Py_BuildValue("d", timestep);
1332}
1333
1334
1335PyObject *compute_fluxes_ext_kinetic(PyObject *self, PyObject *args) {
1336  /*Compute all fluxes and the timestep suitable for all volumes
1337    in domain.
1338
1339    Compute total flux for each conserved quantity using "flux_function_central"
1340
1341    Fluxes across each edge are scaled by edgelengths and summed up
1342    Resulting flux is then scaled by area and stored in
1343    explicit_update for each of the three conserved quantities
1344    stage, xmomentum and ymomentum
1345
1346    The maximal allowable speed computed by the flux_function for each volume
1347    is converted to a timestep that must not be exceeded. The minimum of
1348    those is computed as the next overall timestep.
1349
1350    Python call:
1351    domain.timestep = compute_fluxes(timestep,
1352                                     domain.epsilon,
1353                                     domain.g,
1354                                     domain.neighbours,
1355                                     domain.neighbour_edges,
1356                                     domain.normals,
1357                                     domain.edgelengths,
1358                                     domain.radii,
1359                                     domain.areas,
1360                                     Stage.edge_values,
1361                                     Xmom.edge_values,
1362                                     Ymom.edge_values,
1363                                     Bed.edge_values,
1364                                     Stage.boundary_values,
1365                                     Xmom.boundary_values,
1366                                     Ymom.boundary_values,
1367                                     Stage.explicit_update,
1368                                     Xmom.explicit_update,
1369                                     Ymom.explicit_update,
1370                                     already_computed_flux)
1371
1372
1373    Post conditions:
1374      domain.explicit_update is reset to computed flux values
1375      domain.timestep is set to the largest step satisfying all volumes.
1376
1377
1378  */
1379
1380
1381  PyArrayObject *neighbours, *neighbour_edges,
1382    *normals, *edgelengths, *radii, *areas,
1383    *tri_full_flag,
1384    *stage_edge_values,
1385    *xmom_edge_values,
1386    *ymom_edge_values,
1387    *bed_edge_values,
1388    *stage_boundary_values,
1389    *xmom_boundary_values,
1390    *ymom_boundary_values,
1391    *stage_explicit_update,
1392    *xmom_explicit_update,
1393    *ymom_explicit_update,
1394    *already_computed_flux;//tracks whether the flux across an edge has already been computed
1395
1396
1397  //Local variables
1398  double timestep, max_speed, epsilon, g;
1399  double normal[2], ql[3], qr[3], zl, zr;
1400  double edgeflux[3]; //Work arrays for summing up fluxes
1401
1402  int number_of_elements, k, i, m, n;
1403  int ki, nm=0, ki2; //Index shorthands
1404  static long call=1;
1405
1406
1407  // Convert Python arguments to C
1408  if (!PyArg_ParseTuple(args, "dddOOOOOOOOOOOOOOOOOO",
1409                        &timestep,
1410                        &epsilon,
1411                        &g,
1412                        &neighbours,
1413                        &neighbour_edges,
1414                        &normals,
1415                        &edgelengths, &radii, &areas,
1416                        &tri_full_flag,
1417                        &stage_edge_values,
1418                        &xmom_edge_values,
1419                        &ymom_edge_values,
1420                        &bed_edge_values,
1421                        &stage_boundary_values,
1422                        &xmom_boundary_values,
1423                        &ymom_boundary_values,
1424                        &stage_explicit_update,
1425                        &xmom_explicit_update,
1426                        &ymom_explicit_update,
1427                        &already_computed_flux)) {
1428    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
1429    return NULL;
1430  }
1431  number_of_elements = stage_edge_values -> dimensions[0];
1432  call++;//a static local variable to which already_computed_flux is compared
1433  //set explicit_update to zero for all conserved_quantities.
1434  //This assumes compute_fluxes called before forcing terms
1435  for (k=0; k<number_of_elements; k++) {
1436    ((double *) stage_explicit_update -> data)[k]=0.0;
1437    ((double *) xmom_explicit_update -> data)[k]=0.0;
1438    ((double *) ymom_explicit_update -> data)[k]=0.0;
1439  }
1440  //Loop through neighbours and compute edge flux for each
1441  for (k=0; k<number_of_elements; k++) {
1442    for (i=0; i<3; i++) {
1443      ki = k*3+i;
1444      if (((long *) already_computed_flux->data)[ki]==call)//we've already computed the flux across this edge
1445        continue;
1446      ql[0] = ((double *) stage_edge_values -> data)[ki];
1447      ql[1] = ((double *) xmom_edge_values -> data)[ki];
1448      ql[2] = ((double *) ymom_edge_values -> data)[ki];
1449      zl =    ((double *) bed_edge_values -> data)[ki];
1450
1451      //Quantities at neighbour on nearest face
1452      n = ((long *) neighbours -> data)[ki];
1453      if (n < 0) {
1454        m = -n-1; //Convert negative flag to index
1455        qr[0] = ((double *) stage_boundary_values -> data)[m];
1456        qr[1] = ((double *) xmom_boundary_values -> data)[m];
1457        qr[2] = ((double *) ymom_boundary_values -> data)[m];
1458        zr = zl; //Extend bed elevation to boundary
1459      } else {
1460        m = ((long *) neighbour_edges -> data)[ki];
1461        nm = n*3+m;
1462        qr[0] = ((double *) stage_edge_values -> data)[nm];
1463        qr[1] = ((double *) xmom_edge_values -> data)[nm];
1464        qr[2] = ((double *) ymom_edge_values -> data)[nm];
1465        zr =    ((double *) bed_edge_values -> data)[nm];
1466      }
1467      // Outward pointing normal vector
1468      // normal = domain.normals[k, 2*i:2*i+2]
1469      ki2 = 2*ki; //k*6 + i*2
1470      normal[0] = ((double *) normals -> data)[ki2];
1471      normal[1] = ((double *) normals -> data)[ki2+1];
1472      //Edge flux computation
1473      flux_function_kinetic(ql, qr, zl, zr,
1474                    normal[0], normal[1],
1475                    epsilon, g,
1476                    edgeflux, &max_speed);
1477      //update triangle k
1478      ((long *) already_computed_flux->data)[ki]=call;
1479      ((double *) stage_explicit_update -> data)[k] -= edgeflux[0]*((double *) edgelengths -> data)[ki];
1480      ((double *) xmom_explicit_update -> data)[k] -= edgeflux[1]*((double *) edgelengths -> data)[ki];
1481      ((double *) ymom_explicit_update -> data)[k] -= edgeflux[2]*((double *) edgelengths -> data)[ki];
1482      //update the neighbour n
1483      if (n>=0){
1484        ((long *) already_computed_flux->data)[nm]=call;
1485        ((double *) stage_explicit_update -> data)[n] += edgeflux[0]*((double *) edgelengths -> data)[nm];
1486        ((double *) xmom_explicit_update -> data)[n] += edgeflux[1]*((double *) edgelengths -> data)[nm];
1487        ((double *) ymom_explicit_update -> data)[n] += edgeflux[2]*((double *) edgelengths -> data)[nm];
1488      }
1489      ///for (j=0; j<3; j++) {
1490        ///flux[j] -= edgeflux[j]*((double *) edgelengths -> data)[ki];
1491        ///}
1492        //Update timestep
1493        //timestep = min(timestep, domain.radii[k]/max_speed)
1494        //FIXME: SR Add parameter for CFL condition
1495    if ( ((long *) tri_full_flag -> data)[k] == 1) {
1496            if (max_speed > epsilon) {
1497                timestep = min(timestep, ((double *) radii -> data)[k]/max_speed);
1498                //maxspeed in flux_function is calculated as max(|u+a|,|u-a|)
1499                if (n>=0)
1500                    timestep = min(timestep, ((double *) radii -> data)[n]/max_speed);
1501            }
1502    }
1503    } // end for i
1504    //Normalise by area and store for when all conserved
1505    //quantities get updated
1506    ((double *) stage_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
1507    ((double *) xmom_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
1508    ((double *) ymom_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
1509  } //end for k
1510  return Py_BuildValue("d", timestep);
1511}
1512
1513PyObject *protect(PyObject *self, PyObject *args) {
1514  //
1515  //    protect(minimum_allowed_height, maximum_allowed_speed, wc, zc, xmomc, ymomc)
1516
1517
1518  PyArrayObject
1519  *wc,            //Stage at centroids
1520  *zc,            //Elevation at centroids
1521  *xmomc,         //Momentums at centroids
1522  *ymomc;
1523
1524
1525  int N;
1526  double minimum_allowed_height, maximum_allowed_speed, epsilon;
1527
1528  // Convert Python arguments to C
1529  if (!PyArg_ParseTuple(args, "dddOOOO",
1530                        &minimum_allowed_height,
1531                        &maximum_allowed_speed,
1532                        &epsilon,
1533                        &wc, &zc, &xmomc, &ymomc)) {
1534    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: protect could not parse input arguments");
1535    return NULL;
1536  } 
1537
1538  N = wc -> dimensions[0];
1539
1540  _protect(N,
1541           minimum_allowed_height,
1542           maximum_allowed_speed,
1543           epsilon,
1544           (double*) wc -> data,
1545           (double*) zc -> data,
1546           (double*) xmomc -> data,
1547           (double*) ymomc -> data);
1548
1549  return Py_BuildValue("");
1550}
1551
1552
1553
1554PyObject *balance_deep_and_shallow(PyObject *self, PyObject *args) {
1555  //
1556  //    balance_deep_and_shallow(wc, zc, hc, wv, zv, hv,
1557  //                             xmomc, ymomc, xmomv, ymomv)
1558
1559
1560  PyArrayObject
1561    *wc,            //Stage at centroids
1562    *zc,            //Elevation at centroids
1563    *hc,            //Height at centroids
1564    *wv,            //Stage at vertices
1565    *zv,            //Elevation at vertices
1566    *hv,            //Depths at vertices
1567    *hvbar,         //h-Limited depths at vertices
1568    *xmomc,         //Momentums at centroids and vertices
1569    *ymomc,
1570    *xmomv,
1571    *ymomv;
1572 
1573  PyObject *domain, *Tmp;
1574   
1575  double alpha_balance = 2.0;
1576
1577  int N; //, err;
1578
1579  // Convert Python arguments to C
1580  if (!PyArg_ParseTuple(args, "OOOOOOOOOOOO",
1581                        &domain,
1582                        &wc, &zc, &hc,
1583                        &wv, &zv, &hv, &hvbar,
1584                        &xmomc, &ymomc, &xmomv, &ymomv)) {
1585    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not parse input arguments");
1586    return NULL;
1587  } 
1588         
1589  // Pull out parameters
1590  Tmp = PyObject_GetAttrString(domain, "alpha_balance");
1591  if (!Tmp) {
1592    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not obtain object alpha_balance from domain");
1593    return NULL;
1594  } 
1595  alpha_balance = PyFloat_AsDouble(Tmp);
1596  Py_DECREF(Tmp);
1597
1598
1599  N = wc -> dimensions[0];
1600
1601  _balance_deep_and_shallow(N,
1602                            (double*) wc -> data,
1603                            (double*) zc -> data,
1604                            (double*) hc -> data,
1605                            (double*) wv -> data,
1606                            (double*) zv -> data,
1607                            (double*) hv -> data,
1608                            (double*) hvbar -> data,
1609                                (double*) xmomc -> data,
1610                            (double*) ymomc -> data,
1611                            (double*) xmomv -> data,
1612                            (double*) ymomv -> data,
1613                            alpha_balance);
1614
1615
1616  return Py_BuildValue("");
1617}
1618
1619
1620
1621PyObject *h_limiter(PyObject *self, PyObject *args) {
1622
1623  PyObject *domain, *Tmp;
1624  PyArrayObject
1625    *hv, *hc, //Depth at vertices and centroids
1626    *hvbar,   //Limited depth at vertices (return values)
1627    *neighbours;
1628
1629  int k, i, n, N, k3;
1630  int dimensions[2];
1631  double beta_h; //Safety factor (see config.py)
1632  double *hmin, *hmax, hn;
1633
1634  // Convert Python arguments to C
1635  if (!PyArg_ParseTuple(args, "OOO", &domain, &hc, &hv)) {
1636    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: h_limiter could not parse input arguments");
1637    return NULL;
1638  } 
1639 
1640  neighbours = get_consecutive_array(domain, "neighbours");
1641
1642  //Get safety factor beta_h
1643  Tmp = PyObject_GetAttrString(domain, "beta_h");
1644  if (!Tmp) {
1645    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: h_limiter could not obtain object beta_h from domain");
1646    return NULL;
1647  } 
1648  beta_h = PyFloat_AsDouble(Tmp);
1649
1650  Py_DECREF(Tmp);
1651
1652  N = hc -> dimensions[0];
1653
1654  //Create hvbar
1655  dimensions[0] = N;
1656  dimensions[1] = 3;
1657  hvbar = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
1658
1659
1660  //Find min and max of this and neighbour's centroid values
1661  hmin = malloc(N * sizeof(double));
1662  hmax = malloc(N * sizeof(double));
1663  for (k=0; k<N; k++) {
1664    k3=k*3;
1665
1666    hmin[k] = ((double*) hc -> data)[k];
1667    hmax[k] = hmin[k];
1668
1669    for (i=0; i<3; i++) {
1670      n = ((long*) neighbours -> data)[k3+i];
1671
1672      //Initialise hvbar with values from hv
1673      ((double*) hvbar -> data)[k3+i] = ((double*) hv -> data)[k3+i];
1674
1675      if (n >= 0) {
1676        hn = ((double*) hc -> data)[n]; //Neighbour's centroid value
1677
1678        hmin[k] = min(hmin[k], hn);
1679        hmax[k] = max(hmax[k], hn);
1680      }
1681    }
1682  }
1683
1684  // Call underlying standard routine
1685  _limit(N, beta_h, (double*) hc -> data, (double*) hvbar -> data, hmin, hmax);
1686
1687  // // //Py_DECREF(domain); //FIXME: NEcessary?
1688  free(hmin);
1689  free(hmax);
1690
1691  //return result using PyArray to avoid memory leak
1692  return PyArray_Return(hvbar);
1693  //return Py_BuildValue("");
1694}
1695
1696PyObject *h_limiter_sw(PyObject *self, PyObject *args) {
1697  //a faster version of h_limiter above
1698  PyObject *domain, *Tmp;
1699  PyArrayObject
1700    *hv, *hc, //Depth at vertices and centroids
1701    *hvbar,   //Limited depth at vertices (return values)
1702    *neighbours;
1703
1704  int k, i, N, k3,k0,k1,k2;
1705  int dimensions[2];
1706  double beta_h; //Safety factor (see config.py)
1707  double hmin, hmax, dh[3];
1708// Convert Python arguments to C
1709  if (!PyArg_ParseTuple(args, "OOO", &domain, &hc, &hv)) {
1710    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: h_limiter_sw could not parse input arguments");
1711    return NULL;
1712  } 
1713  neighbours = get_consecutive_array(domain, "neighbours");
1714
1715  //Get safety factor beta_h
1716  Tmp = PyObject_GetAttrString(domain, "beta_h");
1717  if (!Tmp) {
1718    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: h_limiter_sw could not obtain object beta_h from domain");
1719    return NULL;
1720  } 
1721  beta_h = PyFloat_AsDouble(Tmp);
1722
1723  Py_DECREF(Tmp);
1724
1725  N = hc -> dimensions[0];
1726
1727  //Create hvbar
1728  dimensions[0] = N;
1729  dimensions[1] = 3;
1730  hvbar = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
1731  for (k=0;k<N;k++){
1732    k3=k*3;
1733    //get the ids of the neighbours
1734    k0 = ((long*) neighbours -> data)[k3];
1735    k1 = ((long*) neighbours -> data)[k3+1];
1736    k2 = ((long*) neighbours -> data)[k3+2];
1737    //set hvbar provisionally
1738    for (i=0;i<3;i++){
1739      ((double*) hvbar -> data)[k3+i] = ((double*) hv -> data)[k3+i];
1740      dh[i]=((double*) hvbar -> data)[k3+i]-((double*) hc -> data)[k];
1741    }
1742    hmin=((double*) hc -> data)[k];
1743    hmax=hmin;
1744    if (k0>=0){
1745      hmin=min(hmin,((double*) hc -> data)[k0]);
1746      hmax=max(hmax,((double*) hc -> data)[k0]);
1747    }
1748    if (k1>=0){
1749      hmin=min(hmin,((double*) hc -> data)[k1]);
1750      hmax=max(hmax,((double*) hc -> data)[k1]);
1751    }
1752    if (k2>=0){
1753      hmin=min(hmin,((double*) hc -> data)[k2]);
1754      hmax=max(hmax,((double*) hc -> data)[k2]);
1755    }
1756    hmin-=((double*) hc -> data)[k];
1757    hmax-=((double*) hc -> data)[k];
1758    limit_gradient(dh,hmin,hmax,beta_h);
1759    for (i=0;i<3;i++)
1760      ((double*) hvbar -> data)[k3+i] = ((double*) hc -> data)[k]+dh[i];
1761  }
1762  return PyArray_Return(hvbar);
1763}
1764
1765PyObject *assign_windfield_values(PyObject *self, PyObject *args) {
1766  //
1767  //      assign_windfield_values(xmom_update, ymom_update,
1768  //                              s_vec, phi_vec, self.const)
1769
1770
1771
1772  PyArrayObject   //(one element per triangle)
1773  *s_vec,         //Speeds
1774  *phi_vec,       //Bearings
1775  *xmom_update,   //Momentum updates
1776  *ymom_update;
1777
1778
1779  int N;
1780  double cw;
1781
1782  // Convert Python arguments to C
1783  if (!PyArg_ParseTuple(args, "OOOOd",
1784                        &xmom_update,
1785                        &ymom_update,
1786                        &s_vec, &phi_vec,
1787                        &cw)) {
1788    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: assign_windfield_values could not parse input arguments");
1789    return NULL;
1790  } 
1791                       
1792
1793  N = xmom_update -> dimensions[0];
1794
1795  _assign_wind_field_values(N,
1796           (double*) xmom_update -> data,
1797           (double*) ymom_update -> data,
1798           (double*) s_vec -> data,
1799           (double*) phi_vec -> data,
1800           cw);
1801
1802  return Py_BuildValue("");
1803}
1804
1805
1806
1807
1808//////////////////////////////////////////
1809// Method table for python module
1810static struct PyMethodDef MethodTable[] = {
1811  /* The cast of the function is necessary since PyCFunction values
1812   * only take two PyObject* parameters, and rotate() takes
1813   * three.
1814   */
1815
1816  {"rotate", (PyCFunction)rotate, METH_VARARGS | METH_KEYWORDS, "Print out"},
1817  {"extrapolate_second_order_sw", extrapolate_second_order_sw, METH_VARARGS, "Print out"},
1818  {"compute_fluxes_ext_central", compute_fluxes_ext_central, METH_VARARGS, "Print out"},
1819  {"compute_fluxes_ext_kinetic", compute_fluxes_ext_kinetic, METH_VARARGS, "Print out"},
1820  {"gravity", gravity, METH_VARARGS, "Print out"},
1821  {"manning_friction", manning_friction, METH_VARARGS, "Print out"},
1822  {"balance_deep_and_shallow", balance_deep_and_shallow,
1823   METH_VARARGS, "Print out"},
1824  {"h_limiter", h_limiter,
1825   METH_VARARGS, "Print out"},
1826  {"h_limiter_sw", h_limiter_sw,
1827   METH_VARARGS, "Print out"},
1828  {"protect", protect, METH_VARARGS | METH_KEYWORDS, "Print out"},
1829  {"assign_windfield_values", assign_windfield_values,
1830   METH_VARARGS | METH_KEYWORDS, "Print out"},
1831  //{"distribute_to_vertices_and_edges",
1832  // distribute_to_vertices_and_edges, METH_VARARGS},
1833  //{"update_conserved_quantities",
1834  // update_conserved_quantities, METH_VARARGS},
1835  //{"set_initialcondition",
1836  // set_initialcondition, METH_VARARGS},
1837  {NULL, NULL}
1838};
1839
1840// Module initialisation
1841void initshallow_water_ext(void){
1842  Py_InitModule("shallow_water_ext", MethodTable);
1843
1844  import_array();     //Necessary for handling of NumPY structures
1845}
Note: See TracBrowser for help on using the repository browser.