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

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

Implemented optimisation excluding dry cells from flux calculation.
All tests pass and the script run_okushiri_profile.py reports an improvement
in compute_fluxes from 11.25s to 8.58s or 24% faster.
The overall computation was about 40s, so this optimisation improved the
total running time for the problem in question by 7%.

This partially addresses ticket:135

File size: 63.3 KB
Line 
1// Python - C extension module for shallow_water.py
2//
3// To compile (Python2.3):
4//  gcc -c domain_ext.c -I/usr/include/python2.3 -o domain_ext.o -Wall -O
5//  gcc -shared domain_ext.o  -o domain_ext.so
6//
7// or use python compile.py
8//
9// See the module shallow_water.py
10//
11//
12// Ole Nielsen, GA 2004
13
14
15#include "Python.h"
16#include "Numeric/arrayobject.h"
17#include "math.h"
18#include <stdio.h>
19
20//Shared code snippets
21#include "util_ext.h"
22
23const double pi = 3.14159265358979;
24
25// Computational function for rotation
26int _rotate(double *q, double n1, double n2) {
27  /*Rotate the momentum component q (q[1], q[2])
28    from x,y coordinates to coordinates based on normal vector (n1, n2).
29
30    Result is returned in array 3x1 r
31    To rotate in opposite direction, call rotate with (q, n1, -n2)
32
33    Contents of q are changed by this function */
34
35
36  double q1, q2;
37
38  //Shorthands
39  q1 = q[1];  //uh momentum
40  q2 = q[2];  //vh momentum
41
42  //Rotate
43  q[1] =  n1*q1 + n2*q2;
44  q[2] = -n2*q1 + n1*q2;
45
46  return 0;
47}
48
49int find_qmin_and_qmax(double dq0, double dq1, double dq2, double *qmin, double *qmax){
50  //Considering the centroid of an FV triangle and the vertices of its auxiliary triangle, find
51  //qmin=min(q)-qc and qmax=max(q)-qc, where min(q) and max(q) are respectively min and max over the
52  //four values (at the centroid of the FV triangle and the auxiliary triangle vertices),
53  //and qc is the centroid
54  //dq0=q(vertex0)-q(centroid of FV triangle)
55  //dq1=q(vertex1)-q(vertex0)
56  //dq2=q(vertex2)-q(vertex0)
57  if (dq0>=0.0){
58    if (dq1>=dq2){
59      if (dq1>=0.0)
60        *qmax=dq0+dq1;
61      else
62        *qmax=dq0;
63      if ((*qmin=dq0+dq2)<0)
64        ;//qmin is already set to correct value
65      else
66        *qmin=0.0;
67    }
68    else{//dq1<dq2
69      if (dq2>0)
70        *qmax=dq0+dq2;
71      else
72        *qmax=dq0;
73      if ((*qmin=dq0+dq1)<0)
74        ;//qmin is the correct value
75      else
76        *qmin=0.0;
77    }
78  }
79  else{//dq0<0
80    if (dq1<=dq2){
81      if (dq1<0.0)
82        *qmin=dq0+dq1;
83      else
84        *qmin=dq0;
85      if ((*qmax=dq0+dq2)>0.0)
86        ;//qmax is already set to the correct value
87      else
88        *qmax=0.0;
89    }
90    else{//dq1>dq2
91      if (dq2<0.0)
92        *qmin=dq0+dq2;
93      else
94        *qmin=dq0;
95      if ((*qmax=dq0+dq1)>0.0)
96        ;//qmax is already set to the correct value
97      else
98        *qmax=0.0;
99    }
100  }
101  return 0;
102}
103
104int limit_gradient(double *dqv, double qmin, double qmax, double beta_w){
105  //given provisional jumps dqv from the FV triangle centroid to its vertices and
106  //jumps qmin (qmax) between the centroid of the FV triangle and the
107  //minimum (maximum) of the values at the centroid of the FV triangle and the auxiliary triangle vertices,
108  //calculate a multiplicative factor phi by which the provisional vertex jumps are to be limited
109  int i;
110  double r=1000.0, r0=1.0, phi=1.0;
111  static double TINY = 1.0e-100;//to avoid machine accuracy problems.
112  //Any provisional jump with magnitude < TINY does not contribute to the limiting process.
113  for (i=0;i<3;i++){
114    if (dqv[i]<-TINY)
115      r0=qmin/dqv[i];
116    if (dqv[i]>TINY)
117      r0=qmax/dqv[i];
118    r=min(r0,r);
119    //
120  }
121  phi=min(r*beta_w,1.0);
122  for (i=0;i<3;i++)
123    dqv[i]=dqv[i]*phi;
124  return 0;
125}
126
127
128void adjust_froude_number(double *uh,
129                          double h, 
130                          double g) {
131                         
132  // Adjust momentum if Froude number is excessive
133  double max_froude_number = 20.0;                       
134  double froude_number;
135 
136  //Compute Froude number (stability diagnostics)
137  froude_number = *uh/sqrt(g*h)/h;
138
139  if (froude_number > max_froude_number) {
140    printf("---------------------------------------------\n");
141    printf("froude_number=%f (uh=%f, h=%f)\n", froude_number, *uh, h);
142   
143    *uh = *uh/fabs(*uh) * max_froude_number * sqrt(g*h)*h;
144   
145    froude_number = *uh/sqrt(g*h)/h;   
146    printf("Adjusted froude_number=%f (uh=%f, h=%f)\n", froude_number, *uh, h);
147    printf("---------------------------------------------\n");   
148  }
149}
150
151
152
153// Function to obtain speed from momentum and depth.
154// This is used by flux functions
155// Input parameters uh and h may be modified by this function.
156double _compute_speed(double *uh, 
157                      double *h, 
158                      double epsilon, 
159                      double h0) {
160 
161  double u;
162
163  //adjust_froude_number(uh, *h, 9.81); // Highly experimental and
164                                        // probably unneccessary
165 
166  if (*h < epsilon) {
167    *h = 0.0;  //Could have been negative
168    u = 0.0;
169  } else {
170    u = *uh/(*h + h0/ *h);   
171  }
172 
173
174  // Adjust momentum to be consistent with speed
175  *uh = u * *h;
176 
177  return u;
178}
179
180// Computational function for flux computation (using stage w=z+h)
181int flux_function_central(double *q_left, double *q_right,
182                          double z_left, double z_right,
183                          double n1, double n2,
184                          double epsilon, double H0, double g,
185                          double *edgeflux, double *max_speed) {
186
187  /*Compute fluxes between volumes for the shallow water wave equation
188    cast in terms of the 'stage', w = h+z using
189    the 'central scheme' as described in
190
191    Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
192    Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
193    Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.
194
195    The implemented formula is given in equation (3.15) on page 714
196  */
197
198  int i;
199
200  double w_left, h_left, uh_left, vh_left, u_left;
201  double w_right, h_right, uh_right, vh_right, u_right;
202  double v_left, v_right; 
203  double s_min, s_max, soundspeed_left, soundspeed_right;
204  double denom, z;
205  double q_left_copy[3], q_right_copy[3];
206  double flux_right[3], flux_left[3];
207
208  double h0 = H0*H0; //This ensures a good balance when h approaches H0.
209                     //But evidence suggests that h0 can be as little as
210                     //epsilon!
211 
212  //Copy conserved quantities to protect from modification
213  for (i=0; i<3; i++) {
214    q_left_copy[i] = q_left[i];
215    q_right_copy[i] = q_right[i];
216  }
217
218  //Align x- and y-momentum with x-axis
219  _rotate(q_left_copy, n1, n2);
220  _rotate(q_right_copy, n1, n2);
221
222  z = (z_left+z_right)/2; //Take average of field values
223
224  //Compute speeds in x-direction
225  w_left = q_left_copy[0];         
226  h_left = w_left-z;
227  uh_left = q_left_copy[1];
228  u_left = _compute_speed(&uh_left, &h_left, epsilon, h0);
229
230  w_right = q_right_copy[0];
231  h_right = w_right-z;
232  uh_right = q_right_copy[1];
233  u_right = _compute_speed(&uh_right, &h_right, epsilon, h0); 
234
235  //Momentum in y-direction
236  vh_left  = q_left_copy[2];
237  vh_right = q_right_copy[2];
238
239  // Limit y-momentum if necessary 
240  v_left = _compute_speed(&vh_left, &h_left, epsilon, h0);
241  v_right = _compute_speed(&vh_right, &h_right, epsilon, h0);
242
243  //Maximal and minimal wave speeds
244  soundspeed_left  = sqrt(g*h_left);
245  soundspeed_right = sqrt(g*h_right);
246
247  s_max = max(u_left+soundspeed_left, u_right+soundspeed_right);
248  if (s_max < 0.0) s_max = 0.0;
249
250  s_min = min(u_left-soundspeed_left, u_right-soundspeed_right);
251   if (s_min > 0.0) s_min = 0.0;
252
253 
254  //Flux formulas
255  flux_left[0] = u_left*h_left;
256  flux_left[1] = u_left*uh_left + 0.5*g*h_left*h_left;
257  flux_left[2] = u_left*vh_left;
258
259  flux_right[0] = u_right*h_right;
260  flux_right[1] = u_right*uh_right + 0.5*g*h_right*h_right;
261  flux_right[2] = u_right*vh_right;
262
263 
264 
265  //Flux computation
266  denom = s_max-s_min;
267  if (denom == 0.0) {
268    for (i=0; i<3; i++) edgeflux[i] = 0.0;
269    *max_speed = 0.0;
270  } else {
271    for (i=0; i<3; i++) {
272      edgeflux[i] = s_max*flux_left[i] - s_min*flux_right[i];
273      edgeflux[i] += s_max*s_min*(q_right_copy[i]-q_left_copy[i]);
274      edgeflux[i] /= denom;
275    }
276
277    //Maximal wavespeed
278    *max_speed = max(fabs(s_max), fabs(s_min));
279
280    //Rotate back
281    _rotate(edgeflux, n1, -n2);
282  }
283 
284  return 0;
285}
286
287double erfcc(double x){
288    double z,t,result;
289
290    z=fabs(x);
291    t=1.0/(1.0+0.5*z);
292    result=t*exp(-z*z-1.26551223+t*(1.00002368+t*(.37409196+
293         t*(.09678418+t*(-.18628806+t*(.27886807+t*(-1.13520398+
294         t*(1.48851587+t*(-.82215223+t*.17087277)))))))));
295    if (x < 0.0) result = 2.0-result;
296
297    return result;
298    }
299
300
301
302// Computational function for flux computation (using stage w=z+h)
303// FIXME (Ole): Is this used anywhere??
304int flux_function_kinetic(double *q_left, double *q_right,
305                  double z_left, double z_right,
306                  double n1, double n2,
307                  double epsilon, double H0, double g,
308                  double *edgeflux, double *max_speed) {
309
310  /*Compute fluxes between volumes for the shallow water wave equation
311    cast in terms of the 'stage', w = h+z using
312    the 'central scheme' as described in
313
314    Zhang et. al., Advances in Water Resources, 26(6), 2003, 635-647.
315  */
316
317  int i;
318
319  double w_left, h_left, uh_left, vh_left, u_left, F_left;
320  double w_right, h_right, uh_right, vh_right, u_right, F_right;
321  double s_min, s_max, soundspeed_left, soundspeed_right;
322  double z;
323  double q_left_copy[3], q_right_copy[3];
324
325  double h0 = H0*H0; //This ensures a good balance when h approaches H0.
326
327  //Copy conserved quantities to protect from modification
328  for (i=0; i<3; i++) {
329    q_left_copy[i] = q_left[i];
330    q_right_copy[i] = q_right[i];
331  }
332
333  //Align x- and y-momentum with x-axis
334  _rotate(q_left_copy, n1, n2);
335  _rotate(q_right_copy, n1, n2);
336
337  z = (z_left+z_right)/2; //Take average of field values
338
339  //Compute speeds in x-direction
340  w_left = q_left_copy[0];         
341  h_left = w_left-z;
342  uh_left = q_left_copy[1];
343  u_left =_compute_speed(&uh_left, &h_left, epsilon, h0);
344
345  w_right = q_right_copy[0];
346  h_right = w_right-z;
347  uh_right = q_right_copy[1];
348  u_right =_compute_speed(&uh_right, &h_right, epsilon, h0); 
349
350
351  //Momentum in y-direction
352  vh_left  = q_left_copy[2];
353  vh_right = q_right_copy[2];
354
355
356  //Maximal and minimal wave speeds
357  soundspeed_left  = sqrt(g*h_left);
358  soundspeed_right = sqrt(g*h_right);
359
360  s_max = max(u_left+soundspeed_left, u_right+soundspeed_right);
361  if (s_max < 0.0) s_max = 0.0;
362
363  s_min = min(u_left-soundspeed_left, u_right-soundspeed_right);
364  if (s_min > 0.0) s_min = 0.0;
365
366
367  F_left  = 0.0;
368  F_right = 0.0;
369  if (h_left > 0.0) F_left = u_left/sqrt(g*h_left);
370  if (h_right > 0.0) F_right = u_right/sqrt(g*h_right);
371
372  for (i=0; i<3; i++) edgeflux[i] = 0.0;
373  *max_speed = 0.0;
374
375  edgeflux[0] = h_left*u_left/2.0*erfcc(-F_left) +  \
376          h_left*sqrt(g*h_left)/2.0/sqrt(pi)*exp(-(F_left*F_left)) + \
377          h_right*u_right/2.0*erfcc(F_right) -  \
378          h_right*sqrt(g*h_right)/2.0/sqrt(pi)*exp(-(F_right*F_right));
379
380  edgeflux[1] = (h_left*u_left*u_left + g/2.0*h_left*h_left)/2.0*erfcc(-F_left) + \
381          u_left*h_left*sqrt(g*h_left)/2.0/sqrt(pi)*exp(-(F_left*F_left)) + \
382          (h_right*u_right*u_right + g/2.0*h_right*h_right)/2.0*erfcc(F_right) -  \
383          u_right*h_right*sqrt(g*h_right)/2.0/sqrt(pi)*exp(-(F_right*F_right));
384
385  edgeflux[2] = vh_left*u_left/2.0*erfcc(-F_left) + \
386          vh_left*sqrt(g*h_left)/2.0/sqrt(pi)*exp(-(F_left*F_left)) + \
387          vh_right*u_right/2.0*erfcc(F_right) - \
388          vh_right*sqrt(g*h_right)/2.0/sqrt(pi)*exp(-(F_right*F_right));
389
390  //Maximal wavespeed
391  *max_speed = max(fabs(s_max), fabs(s_min));
392
393  //Rotate back
394  _rotate(edgeflux, n1, -n2);
395
396  return 0;
397}
398
399
400
401
402void _manning_friction(double g, double eps, int N,
403                       double* w, double* z,
404                       double* uh, double* vh,
405                       double* eta, double* xmom, double* ymom) {
406
407  int k;
408  double S, h;
409
410  for (k=0; k<N; k++) {
411    if (eta[k] > eps) {
412      h = w[k]-z[k];
413      if (h >= eps) {
414        S = -g * eta[k]*eta[k] * sqrt((uh[k]*uh[k] + vh[k]*vh[k]));
415        S /= pow(h, 7.0/3);      //Expensive (on Ole's home computer)
416        //S /= exp(7.0/3.0*log(h));      //seems to save about 15% over manning_friction
417        //S /= h*h*(1 + h/3.0 - h*h/9.0); //FIXME: Could use a Taylor expansion
418
419
420        //Update momentum
421        xmom[k] += S*uh[k];
422        ymom[k] += S*vh[k];
423      }
424    }
425  }
426}
427
428
429/*
430void _manning_friction_explicit(double g, double eps, int N,
431                       double* w, double* z,
432                       double* uh, double* vh,
433                       double* eta, double* xmom, double* ymom) {
434
435  int k;
436  double S, h;
437
438  for (k=0; k<N; k++) {
439    if (eta[k] > eps) {
440      h = w[k]-z[k];
441      if (h >= eps) {
442        S = -g * eta[k]*eta[k] * sqrt((uh[k]*uh[k] + vh[k]*vh[k]));
443        S /= pow(h, 7.0/3);      //Expensive (on Ole's home computer)
444        //S /= exp(7.0/3.0*log(h));      //seems to save about 15% over manning_friction
445        //S /= h*h*(1 + h/3.0 - h*h/9.0); //FIXME: Could use a Taylor expansion
446
447
448        //Update momentum
449        xmom[k] += S*uh[k];
450        ymom[k] += S*vh[k];
451      }
452    }
453  }
454}
455*/
456
457int _balance_deep_and_shallow(int N,
458                              double* wc,
459                              double* zc,
460                              double* hc,
461                              double* wv,
462                              double* zv,
463                              double* hv,
464                              double* hvbar,
465                              double* xmomc,
466                              double* ymomc,
467                              double* xmomv,
468                              double* ymomv,
469                              double H0,
470                              int tight_slope_limiters,
471                              double alpha_balance) {
472
473  int k, k3, i;
474  double dz, hmin, alpha, h_diff;
475
476  //Compute linear combination between w-limited stages and
477  //h-limited stages close to the bed elevation.
478
479  for (k=0; k<N; k++) {
480    // Compute maximal variation in bed elevation
481    // This quantitiy is
482    //     dz = max_i abs(z_i - z_c)
483    // and it is independent of dimension
484    // In the 1d case zc = (z0+z1)/2
485    // In the 2d case zc = (z0+z1+z2)/3
486
487    k3 = 3*k;
488
489    //FIXME: Try with this one precomputed
490    dz = 0.0;
491    hmin = hv[k3];
492    for (i=0; i<3; i++) {
493      if (tight_slope_limiters == 0) { 
494        dz = max(dz, fabs(zv[k3+i]-zc[k]));
495      }
496     
497      hmin = min(hmin, hv[k3+i]);
498    }
499
500
501    //Create alpha in [0,1], where alpha==0 means using the h-limited
502    //stage and alpha==1 means using the w-limited stage as
503    //computed by the gradient limiter (both 1st or 2nd order)
504   
505   
506    if (tight_slope_limiters == 0) {     
507      //If hmin > dz/alpha_balance then alpha = 1 and the bed will have no
508      //effect
509      //If hmin < 0 then alpha = 0 reverting to constant height above bed.
510      //The parameter alpha_balance==2 by default
511
512     
513      if (dz > 0.0) {
514        alpha = max( min( alpha_balance*hmin/dz, 1.0), 0.0 );     
515      } else {
516        alpha = 1.0;  //Flat bed
517      }
518      //printf("Using old style limiter\n");
519     
520    } else {
521
522      // 2007 Balanced Limiter
523   
524      // Make alpha as large as possible but still ensure that
525      // final depth is positive
526   
527      if (hmin < H0) {
528        alpha = 1.0;
529        for (i=0; i<3; i++) {
530
531          h_diff = hvbar[k3+i] - hv[k3+i];
532       
533          if (h_diff <= 0) {
534            // Deep water triangle is further away from bed than
535            // shallow water (hbar < h). Any alpha will do
536         
537          } else { 
538            // Denominator is positive which means that we need some of the
539            // h-limited stage.
540           
541            alpha = min(alpha, (hvbar[k3+i] - H0)/h_diff);
542          }
543        }
544
545        // Ensure alpha in [0,1]
546        if (alpha>1.0) alpha=1.0;
547        if (alpha<0.0) alpha=0.0;
548       
549      } else {
550        // Use w-limited stage exclusively
551        alpha = 1.0;       
552      }
553    }
554       
555   
556       
557    //printf("k=%d, hmin=%.2f, dz=%.2f, alpha=%.2f, alpha_balance=%.2f\n",
558    //     k, hmin, dz, alpha, alpha_balance);
559
560    //printf("dz = %.3f, alpha = %.8f\n", dz, alpha);
561
562    //  Let
563    //
564    //    wvi be the w-limited stage (wvi = zvi + hvi)
565    //    wvi- be the h-limited state (wvi- = zvi + hvi-)
566    //
567    //
568    //  where i=0,1,2 denotes the vertex ids
569    //
570    //  Weighted balance between w-limited and h-limited stage is
571    //
572    //    wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi)
573    //
574    //  It follows that the updated wvi is
575    //    wvi := zvi + (1-alpha)*hvi- + alpha*hvi
576    //
577    //   Momentum is balanced between constant and limited
578
579    if (alpha < 1) {
580      for (i=0; i<3; i++) {
581        wv[k3+i] = zv[k3+i] + (1-alpha)*hvbar[k3+i] + alpha*hv[k3+i];
582
583        //Update momentum as a linear combination of
584        //xmomc and ymomc (shallow) and momentum
585        //from extrapolator xmomv and ymomv (deep).
586        //FIXME (Ole): Is this really needed?
587        xmomv[k3+i] = (1-alpha)*xmomc[k] + alpha*xmomv[k3+i];
588        ymomv[k3+i] = (1-alpha)*ymomc[k] + alpha*ymomv[k3+i];
589      }
590    }
591  }
592  return 0;
593}
594
595
596
597int _protect(int N,
598             double minimum_allowed_height,
599             double maximum_allowed_speed,
600             double epsilon,
601             double* wc,
602             double* zc,
603             double* xmomc,
604             double* ymomc) {
605
606  int k;
607  double hc;
608  double u, v, reduced_speed;
609
610  //Protect against initesimal and negative heights
611  for (k=0; k<N; k++) {
612    hc = wc[k] - zc[k];
613
614    if (hc < minimum_allowed_height) {
615       
616      //Old code: Set momentum to zero and ensure h is non negative
617      //xmomc[k] = 0.0;
618      //ymomc[k] = 0.0;
619      //if (hc <= 0.0) wc[k] = zc[k];
620
621
622      //New code: Adjust momentum to guarantee speeds are physical
623      //          ensure h is non negative
624      //FIXME (Ole): This is only implemented in this C extension and
625      //             has no Python equivalent
626           
627      if (hc <= 0.0) {
628        wc[k] = zc[k];
629        xmomc[k] = 0.0;
630        ymomc[k] = 0.0;
631      } else {
632        //Reduce excessive speeds derived from division by small hc
633        //FIXME (Ole): This may be unnecessary with new slope limiters
634        //in effect.
635       
636        u = xmomc[k]/hc;
637        if (fabs(u) > maximum_allowed_speed) {
638          reduced_speed = maximum_allowed_speed * u/fabs(u);
639          //printf("Speed (u) has been reduced from %.3f to %.3f\n",
640          //     u, reduced_speed);
641          xmomc[k] = reduced_speed * hc;
642        }
643
644        v = ymomc[k]/hc;
645        if (fabs(v) > maximum_allowed_speed) {
646          reduced_speed = maximum_allowed_speed * v/fabs(v);
647          //printf("Speed (v) has been reduced from %.3f to %.3f\n",
648          //     v, reduced_speed);
649          ymomc[k] = reduced_speed * hc;
650        }
651      }
652    }
653  }
654  return 0;
655}
656
657
658
659
660int _assign_wind_field_values(int N,
661                              double* xmom_update,
662                              double* ymom_update,
663                              double* s_vec,
664                              double* phi_vec,
665                              double cw) {
666
667  //Assign windfield values to momentum updates
668
669  int k;
670  double S, s, phi, u, v;
671
672  for (k=0; k<N; k++) {
673
674    s = s_vec[k];
675    phi = phi_vec[k];
676
677    //Convert to radians
678    phi = phi*pi/180;
679
680    //Compute velocity vector (u, v)
681    u = s*cos(phi);
682    v = s*sin(phi);
683
684    //Compute wind stress
685    S = cw * sqrt(u*u + v*v);
686    xmom_update[k] += S*u;
687    ymom_update[k] += S*v;
688  }
689  return 0;
690}
691
692
693
694///////////////////////////////////////////////////////////////////
695// Gateways to Python
696
697PyObject *gravity(PyObject *self, PyObject *args) {
698  //
699  //  gravity(g, h, v, x, xmom, ymom)
700  //
701
702
703  PyArrayObject *h, *v, *x, *xmom, *ymom;
704  int k, i, N, k3, k6;
705  double g, avg_h, zx, zy;
706  double x0, y0, x1, y1, x2, y2, z0, z1, z2;
707
708  if (!PyArg_ParseTuple(args, "dOOOOO",
709                        &g, &h, &v, &x,
710                        &xmom, &ymom)) {
711    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: gravity could not parse input arguments");
712    return NULL;
713  }
714
715  N = h -> dimensions[0];
716  for (k=0; k<N; k++) {
717    k3 = 3*k;  // base index
718    k6 = 6*k;  // base index
719
720    avg_h = 0.0;
721    for (i=0; i<3; i++) {
722      avg_h += ((double *) h -> data)[k3+i];
723    }
724    avg_h /= 3;
725
726
727    //Compute bed slope
728    x0 = ((double*) x -> data)[k6 + 0];
729    y0 = ((double*) x -> data)[k6 + 1];
730    x1 = ((double*) x -> data)[k6 + 2];
731    y1 = ((double*) x -> data)[k6 + 3];
732    x2 = ((double*) x -> data)[k6 + 4];
733    y2 = ((double*) x -> data)[k6 + 5];
734
735
736    z0 = ((double*) v -> data)[k3 + 0];
737    z1 = ((double*) v -> data)[k3 + 1];
738    z2 = ((double*) v -> data)[k3 + 2];
739
740    _gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2, &zx, &zy);
741
742    //Update momentum
743    ((double*) xmom -> data)[k] += -g*zx*avg_h;
744    ((double*) ymom -> data)[k] += -g*zy*avg_h;
745  }
746
747  return Py_BuildValue("");
748}
749
750
751PyObject *manning_friction(PyObject *self, PyObject *args) {
752  //
753  // manning_friction(g, eps, h, uh, vh, eta, xmom_update, ymom_update)
754  //
755
756
757  PyArrayObject *w, *z, *uh, *vh, *eta, *xmom, *ymom;
758  int N;
759  double g, eps;
760
761  if (!PyArg_ParseTuple(args, "ddOOOOOOO",
762                        &g, &eps, &w, &z, &uh, &vh, &eta,
763                        &xmom, &ymom)) {
764    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: manning_friction could not parse input arguments");
765    return NULL;
766  }
767
768
769  N = w -> dimensions[0];
770  _manning_friction(g, eps, N,
771                    (double*) w -> data,
772                    (double*) z -> data,
773                    (double*) uh -> data,
774                    (double*) vh -> data,
775                    (double*) eta -> data,
776                    (double*) xmom -> data,
777                    (double*) ymom -> data);
778
779  return Py_BuildValue("");
780}
781
782
783/*
784PyObject *manning_friction_explicit(PyObject *self, PyObject *args) {
785  //
786  // manning_friction_explicit(g, eps, h, uh, vh, eta, xmom_update, ymom_update)
787  //
788
789
790  PyArrayObject *w, *z, *uh, *vh, *eta, *xmom, *ymom;
791  int N;
792  double g, eps;
793
794  if (!PyArg_ParseTuple(args, "ddOOOOOOO",
795                        &g, &eps, &w, &z, &uh, &vh, &eta,
796                        &xmom, &ymom))
797    return NULL;
798
799  N = w -> dimensions[0];
800  _manning_friction_explicit(g, eps, N,
801                    (double*) w -> data,
802                    (double*) z -> data,
803                    (double*) uh -> data,
804                    (double*) vh -> data,
805                    (double*) eta -> data,
806                    (double*) xmom -> data,
807                    (double*) ymom -> data);
808
809  return Py_BuildValue("");
810}
811*/
812
813PyObject *extrapolate_second_order_sw(PyObject *self, PyObject *args) {
814  /*Compute the vertex values based on a linear reconstruction on each triangle
815    These values are calculated as follows:
816    1) For each triangle not adjacent to a boundary, we consider the auxiliary triangle
817    formed by the centroids of its three neighbours.
818    2) For each conserved quantity, we integrate around the auxiliary triangle's boundary the product
819    of the quantity and the outward normal vector. Dividing by the triangle area gives (a,b), the average
820    of the vector (q_x,q_y) on the auxiliary triangle. We suppose that the linear reconstruction on the
821    original triangle has gradient (a,b).
822    3) Provisional vertex jumps dqv[0,1,2] are computed and these are then limited by calling the functions
823    find_qmin_and_qmax and limit_gradient
824
825    Python call:
826    extrapolate_second_order_sw(domain.surrogate_neighbours,
827                                domain.number_of_boundaries
828                                domain.centroid_coordinates,
829                                Stage.centroid_values
830                                Xmom.centroid_values
831                                Ymom.centroid_values
832                                domain.vertex_coordinates,
833                                Stage.vertex_values,
834                                Xmom.vertex_values,
835                                Ymom.vertex_values)
836
837    Post conditions:
838            The vertices of each triangle have values from a limited linear reconstruction
839            based on centroid values
840
841  */
842  PyArrayObject *surrogate_neighbours,
843    *number_of_boundaries,
844    *centroid_coordinates,
845    *stage_centroid_values,
846    *xmom_centroid_values,
847    *ymom_centroid_values,
848        *elevation_centroid_values,
849    *vertex_coordinates,
850    *stage_vertex_values,
851    *xmom_vertex_values,
852    *ymom_vertex_values,
853        *elevation_vertex_values;
854  PyObject *domain, *Tmp;
855  //Local variables
856  double a, b;//gradient vector, not stored but used to calculate vertex values from centroids
857  int number_of_elements,k,k0,k1,k2,k3,k6,coord_index,i;
858  double x,y,x0,y0,x1,y1,x2,y2,xv0,yv0,xv1,yv1,xv2,yv2;//vertices of the auxiliary triangle
859  double dx1,dx2,dy1,dy2,dxv0,dxv1,dxv2,dyv0,dyv1,dyv2,dq0,dq1,dq2,area2;
860  double dqv[3], qmin, qmax, hmin;
861  double hc, h0, h1, h2;
862  double beta_w, beta_w_dry, beta_uh, beta_uh_dry, beta_vh, beta_vh_dry, beta_tmp;
863  double minimum_allowed_height;
864  //provisional jumps from centroids to v'tices and safety factor re limiting
865  //by which these jumps are limited
866  // Convert Python arguments to C
867  if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOO",
868                        &domain,
869                        &surrogate_neighbours,
870                        &number_of_boundaries,
871                        &centroid_coordinates,
872                        &stage_centroid_values,
873                        &xmom_centroid_values,
874                        &ymom_centroid_values,
875                        &elevation_centroid_values,
876                        &vertex_coordinates,
877                        &stage_vertex_values,
878                        &xmom_vertex_values,
879                        &ymom_vertex_values,
880                        &elevation_vertex_values)) {
881    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
882    return NULL;
883  }
884
885  //get the safety factor beta_w, set in the config.py file. This is used in the limiting process
886  Tmp = PyObject_GetAttrString(domain, "beta_w");
887  if (!Tmp) {
888    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_w from domain");
889    return NULL;
890  } 
891  beta_w = PyFloat_AsDouble(Tmp);
892  Py_DECREF(Tmp);
893 
894  Tmp = PyObject_GetAttrString(domain, "beta_w_dry");
895  if (!Tmp) {
896    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_w_dry from domain");
897    return NULL;
898  } 
899  beta_w_dry = PyFloat_AsDouble(Tmp);
900  Py_DECREF(Tmp);
901 
902  Tmp = PyObject_GetAttrString(domain, "beta_uh");
903  if (!Tmp) {
904    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_uh from domain");
905    return NULL;
906  } 
907  beta_uh = PyFloat_AsDouble(Tmp);
908  Py_DECREF(Tmp);
909 
910  Tmp = PyObject_GetAttrString(domain, "beta_uh_dry");
911  if (!Tmp) {
912    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_uh_dry from domain");
913    return NULL;
914  } 
915  beta_uh_dry = PyFloat_AsDouble(Tmp);
916  Py_DECREF(Tmp); 
917
918  Tmp = PyObject_GetAttrString(domain, "beta_vh");
919  if (!Tmp) {
920    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_vh from domain");
921    return NULL;
922  } 
923  beta_vh = PyFloat_AsDouble(Tmp);
924  Py_DECREF(Tmp);
925 
926  Tmp = PyObject_GetAttrString(domain, "beta_vh_dry");
927  if (!Tmp) {
928    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_vh_dry from domain");
929    return NULL;
930  } 
931  beta_vh_dry = PyFloat_AsDouble(Tmp);
932  Py_DECREF(Tmp);
933 
934  Tmp = PyObject_GetAttrString(domain, "minimum_allowed_height");
935  if (!Tmp) {
936    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object minimum_allowed_heigt");
937    return NULL;
938  } 
939  minimum_allowed_height = PyFloat_AsDouble(Tmp);
940  Py_DECREF(Tmp); 
941 
942  number_of_elements = stage_centroid_values -> dimensions[0];
943  for (k=0; k<number_of_elements; k++) {
944    k3=k*3;
945    k6=k*6;
946
947    if (((long *) number_of_boundaries->data)[k]==3){/*no neighbours, set gradient on the triangle to zero*/
948      ((double *) stage_vertex_values->data)[k3]=((double *)stage_centroid_values->data)[k];
949      ((double *) stage_vertex_values->data)[k3+1]=((double *)stage_centroid_values->data)[k];
950      ((double *) stage_vertex_values->data)[k3+2]=((double *)stage_centroid_values->data)[k];
951      ((double *) xmom_vertex_values->data)[k3]=((double *)xmom_centroid_values->data)[k];
952      ((double *) xmom_vertex_values->data)[k3+1]=((double *)xmom_centroid_values->data)[k];
953      ((double *) xmom_vertex_values->data)[k3+2]=((double *)xmom_centroid_values->data)[k];
954      ((double *) ymom_vertex_values->data)[k3]=((double *)ymom_centroid_values->data)[k];
955      ((double *) ymom_vertex_values->data)[k3+1]=((double *)ymom_centroid_values->data)[k];
956      ((double *) ymom_vertex_values->data)[k3+2]=((double *)ymom_centroid_values->data)[k];
957      continue;
958    }
959    else {//we will need centroid coordinates and vertex coordinates of the triangle
960      //get the vertex coordinates of the FV triangle
961      xv0=((double *)vertex_coordinates->data)[k6]; yv0=((double *)vertex_coordinates->data)[k6+1];
962      xv1=((double *)vertex_coordinates->data)[k6+2]; yv1=((double *)vertex_coordinates->data)[k6+3];
963      xv2=((double *)vertex_coordinates->data)[k6+4]; yv2=((double *)vertex_coordinates->data)[k6+5];
964      //get the centroid coordinates of the FV triangle
965      coord_index=2*k;
966      x=((double *)centroid_coordinates->data)[coord_index];
967      y=((double *)centroid_coordinates->data)[coord_index+1];
968      //store x- and y- differentials for the vertices of the FV triangle relative to the centroid
969      dxv0=xv0-x; dxv1=xv1-x; dxv2=xv2-x;
970      dyv0=yv0-y; dyv1=yv1-y; dyv2=yv2-y;
971    }
972    if (((long *)number_of_boundaries->data)[k]<=1){
973      //if no boundaries, auxiliary triangle is formed from the centroids of the three neighbours
974      //if one boundary, auxiliary triangle is formed from this centroid and its two neighbours
975      k0=((long *)surrogate_neighbours->data)[k3];
976      k1=((long *)surrogate_neighbours->data)[k3+1];
977      k2=((long *)surrogate_neighbours->data)[k3+2];
978      //get the auxiliary triangle's vertex coordinates (really the centroids of neighbouring triangles)
979      coord_index=2*k0;
980      x0=((double *)centroid_coordinates->data)[coord_index];
981      y0=((double *)centroid_coordinates->data)[coord_index+1];
982      coord_index=2*k1;
983      x1=((double *)centroid_coordinates->data)[coord_index];
984      y1=((double *)centroid_coordinates->data)[coord_index+1];
985      coord_index=2*k2;
986      x2=((double *)centroid_coordinates->data)[coord_index];
987      y2=((double *)centroid_coordinates->data)[coord_index+1];
988      //store x- and y- differentials for the vertices of the auxiliary triangle
989      dx1=x1-x0; dx2=x2-x0;
990      dy1=y1-y0; dy2=y2-y0;
991      //calculate 2*area of the auxiliary triangle
992      area2 = dy2*dx1 - dy1*dx2;//the triangle is guaranteed to be counter-clockwise
993      //If the mesh is 'weird' near the boundary, the trianlge might be flat or clockwise:
994      if (area2<=0) {
995        PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: negative triangle area encountered");
996        return NULL;
997      } 
998     
999
1000      //### Calculate heights of neighbouring cells
1001      hc = ((double *)stage_centroid_values->data)[k]  - ((double *)elevation_centroid_values->data)[k];
1002      h0 = ((double *)stage_centroid_values->data)[k0] - ((double *)elevation_centroid_values->data)[k0];
1003      h1 = ((double *)stage_centroid_values->data)[k1] - ((double *)elevation_centroid_values->data)[k1];
1004      h2 = ((double *)stage_centroid_values->data)[k2] - ((double *)elevation_centroid_values->data)[k2];
1005      hmin = min(hc,min(h0,min(h1,h2)));
1006     
1007      //### stage ###
1008      //calculate the difference between vertex 0 of the auxiliary triangle and the FV triangle centroid
1009      dq0=((double *)stage_centroid_values->data)[k0]-((double *)stage_centroid_values->data)[k];
1010      //calculate differentials between the vertices of the auxiliary triangle
1011      dq1=((double *)stage_centroid_values->data)[k1]-((double *)stage_centroid_values->data)[k0];
1012      dq2=((double *)stage_centroid_values->data)[k2]-((double *)stage_centroid_values->data)[k0];
1013      //calculate the gradient of stage on the auxiliary triangle
1014      a = dy2*dq1 - dy1*dq2;
1015      a /= area2;
1016      b = dx1*dq2 - dx2*dq1;
1017      b /= area2;
1018      //calculate provisional jumps in stage from the centroid of the FV tri to its vertices, to be limited
1019      dqv[0]=a*dxv0+b*dyv0;
1020      dqv[1]=a*dxv1+b*dyv1;
1021      dqv[2]=a*dxv2+b*dyv2;
1022      //now we want to find min and max of the centroid and the vertices of the auxiliary triangle
1023      //and compute jumps from the centroid to the min and max
1024      find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
1025      // Playing with dry wet interface
1026      hmin = qmin;
1027      beta_tmp = beta_w;
1028      if (hmin<minimum_allowed_height)
1029        beta_tmp = beta_w_dry;
1030      limit_gradient(dqv,qmin,qmax,beta_tmp);//the gradient will be limited
1031      for (i=0;i<3;i++)
1032        ((double *)stage_vertex_values->data)[k3+i]=((double *)stage_centroid_values->data)[k]+dqv[i];
1033     
1034      //### xmom ###
1035      //calculate the difference between vertex 0 of the auxiliary triangle and the FV triangle centroid
1036      dq0=((double *)xmom_centroid_values->data)[k0]-((double *)xmom_centroid_values->data)[k];
1037      //calculate differentials between the vertices of the auxiliary triangle
1038      dq1=((double *)xmom_centroid_values->data)[k1]-((double *)xmom_centroid_values->data)[k0];
1039      dq2=((double *)xmom_centroid_values->data)[k2]-((double *)xmom_centroid_values->data)[k0];
1040      //calculate the gradient of xmom on the auxiliary triangle
1041      a = dy2*dq1 - dy1*dq2;
1042      a /= area2;
1043      b = dx1*dq2 - dx2*dq1;
1044      b /= area2;
1045      //calculate provisional jumps in stage from the centroid of the FV tri to its vertices, to be limited
1046      dqv[0]=a*dxv0+b*dyv0;
1047      dqv[1]=a*dxv1+b*dyv1;
1048      dqv[2]=a*dxv2+b*dyv2;
1049      //now we want to find min and max of the centroid and the vertices of the auxiliary triangle
1050      //and compute jumps from the centroid to the min and max
1051      find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
1052      beta_tmp = beta_uh;
1053      if (hmin<minimum_allowed_height)
1054        beta_tmp = beta_uh_dry;
1055      limit_gradient(dqv,qmin,qmax,beta_tmp);//the gradient will be limited
1056      for (i=0;i<3;i++)
1057        ((double *)xmom_vertex_values->data)[k3+i]=((double *)xmom_centroid_values->data)[k]+dqv[i];
1058     
1059      //### ymom ###
1060      //calculate the difference between vertex 0 of the auxiliary triangle and the FV triangle centroid
1061      dq0=((double *)ymom_centroid_values->data)[k0]-((double *)ymom_centroid_values->data)[k];
1062      //calculate differentials between the vertices of the auxiliary triangle
1063      dq1=((double *)ymom_centroid_values->data)[k1]-((double *)ymom_centroid_values->data)[k0];
1064      dq2=((double *)ymom_centroid_values->data)[k2]-((double *)ymom_centroid_values->data)[k0];
1065      //calculate the gradient of xmom on the auxiliary triangle
1066      a = dy2*dq1 - dy1*dq2;
1067      a /= area2;
1068      b = dx1*dq2 - dx2*dq1;
1069      b /= area2;
1070      //calculate provisional jumps in stage from the centroid of the FV tri to its vertices, to be limited
1071      dqv[0]=a*dxv0+b*dyv0;
1072      dqv[1]=a*dxv1+b*dyv1;
1073      dqv[2]=a*dxv2+b*dyv2;
1074      //now we want to find min and max of the centroid and the vertices of the auxiliary triangle
1075      //and compute jumps from the centroid to the min and max
1076      find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
1077      beta_tmp = beta_vh;
1078      if (hmin<minimum_allowed_height)
1079        beta_tmp = beta_vh_dry;
1080      limit_gradient(dqv,qmin,qmax,beta_tmp);//the gradient will be limited
1081      for (i=0;i<3;i++)
1082        ((double *)ymom_vertex_values->data)[k3+i]=((double *)ymom_centroid_values->data)[k]+dqv[i];
1083    }//if (number_of_boundaries[k]<=1)
1084    else{//number_of_boundaries==2
1085      //one internal neighbour and gradient is in direction of the neighbour's centroid
1086      //find the only internal neighbour
1087      for (k2=k3;k2<k3+3;k2++){//k2 just indexes the edges of triangle k
1088        if (((long *)surrogate_neighbours->data)[k2]!=k)//find internal neighbour of triabngle k
1089          break;
1090      }
1091      if ((k2==k3+3)) {//if we didn't find an internal neighbour
1092        PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: Internal neighbour not found");     
1093        return NULL;//error
1094      }
1095     
1096      k1=((long *)surrogate_neighbours->data)[k2];
1097      //the coordinates of the triangle are already (x,y). Get centroid of the neighbour (x1,y1)
1098      coord_index=2*k1;
1099      x1=((double *)centroid_coordinates->data)[coord_index];
1100      y1=((double *)centroid_coordinates->data)[coord_index+1];
1101      //compute x- and y- distances between the centroid of the FV triangle and that of its neighbour
1102      dx1=x1-x; dy1=y1-y;
1103      //set area2 as the square of the distance
1104      area2=dx1*dx1+dy1*dy1;
1105      //set dx2=(x1-x0)/((x1-x0)^2+(y1-y0)^2) and dy2=(y1-y0)/((x1-x0)^2+(y1-y0)^2) which
1106      //respectively correspond to the x- and y- gradients of the conserved quantities
1107      dx2=1.0/area2;
1108      dy2=dx2*dy1;
1109      dx2*=dx1;
1110     
1111      //## stage ###
1112      //compute differentials
1113      dq1=((double *)stage_centroid_values->data)[k1]-((double *)stage_centroid_values->data)[k];
1114      //calculate the gradient between the centroid of the FV triangle and that of its neighbour
1115      a=dq1*dx2;
1116      b=dq1*dy2;
1117      //calculate provisional vertex jumps, to be limited
1118      dqv[0]=a*dxv0+b*dyv0;
1119      dqv[1]=a*dxv1+b*dyv1;
1120      dqv[2]=a*dxv2+b*dyv2;
1121      //now limit the jumps
1122      if (dq1>=0.0){
1123        qmin=0.0;
1124        qmax=dq1;
1125      }
1126      else{
1127        qmin=dq1;
1128        qmax=0.0;
1129      }
1130     
1131     
1132      limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
1133      for (i=0;i<3;i++)
1134        ((double *)stage_vertex_values->data)[k3+i]=((double *)stage_centroid_values->data)[k]+dqv[i];
1135     
1136      //## xmom ###
1137      //compute differentials
1138      dq1=((double *)xmom_centroid_values->data)[k1]-((double *)xmom_centroid_values->data)[k];
1139      //calculate the gradient between the centroid of the FV triangle and that of its neighbour
1140      a=dq1*dx2;
1141      b=dq1*dy2;
1142      //calculate provisional vertex jumps, to be limited
1143      dqv[0]=a*dxv0+b*dyv0;
1144      dqv[1]=a*dxv1+b*dyv1;
1145      dqv[2]=a*dxv2+b*dyv2;
1146      //now limit the jumps
1147      if (dq1>=0.0){
1148        qmin=0.0;
1149        qmax=dq1;
1150      }
1151      else{
1152        qmin=dq1;
1153        qmax=0.0;
1154      }
1155      limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
1156      for (i=0;i<3;i++)
1157        ((double *)xmom_vertex_values->data)[k3+i]=((double *)xmom_centroid_values->data)[k]+dqv[i];
1158     
1159      //## ymom ###
1160      //compute differentials
1161      dq1=((double *)ymom_centroid_values->data)[k1]-((double *)ymom_centroid_values->data)[k];
1162      //calculate the gradient between the centroid of the FV triangle and that of its neighbour
1163      a=dq1*dx2;
1164      b=dq1*dy2;
1165      //calculate provisional vertex jumps, to be limited
1166      dqv[0]=a*dxv0+b*dyv0;
1167      dqv[1]=a*dxv1+b*dyv1;
1168      dqv[2]=a*dxv2+b*dyv2;
1169      //now limit the jumps
1170      if (dq1>=0.0){
1171        qmin=0.0;
1172        qmax=dq1;
1173      }
1174      else{
1175        qmin=dq1;
1176        qmax=0.0;
1177      }
1178      limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
1179      for (i=0;i<3;i++)
1180        ((double *)ymom_vertex_values->data)[k3+i]=((double *)ymom_centroid_values->data)[k]+dqv[i];
1181    }//else [number_of_boudaries==2]
1182  }//for k=0 to number_of_elements-1
1183  return Py_BuildValue("");
1184}//extrapolate_second-order_sw
1185
1186
1187PyObject *rotate(PyObject *self, PyObject *args, PyObject *kwargs) {
1188  //
1189  // r = rotate(q, normal, direction=1)
1190  //
1191  // Where q is assumed to be a Float numeric array of length 3 and
1192  // normal a Float numeric array of length 2.
1193
1194  // FIXME(Ole): I don't think this is used anymore
1195
1196  PyObject *Q, *Normal;
1197  PyArrayObject *q, *r, *normal;
1198
1199  static char *argnames[] = {"q", "normal", "direction", NULL};
1200  int dimensions[1], i, direction=1;
1201  double n1, n2;
1202
1203  // Convert Python arguments to C
1204  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OO|i", argnames,
1205                                   &Q, &Normal, &direction)) {
1206    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: rotate could not parse input arguments");
1207    return NULL;
1208  } 
1209
1210  //Input checks (convert sequences into numeric arrays)
1211  q = (PyArrayObject *)
1212    PyArray_ContiguousFromObject(Q, PyArray_DOUBLE, 0, 0);
1213  normal = (PyArrayObject *)
1214    PyArray_ContiguousFromObject(Normal, PyArray_DOUBLE, 0, 0);
1215
1216
1217  if (normal -> dimensions[0] != 2) {
1218    PyErr_SetString(PyExc_RuntimeError, "Normal vector must have 2 components");
1219    return NULL;
1220  }
1221
1222  //Allocate space for return vector r (don't DECREF)
1223  dimensions[0] = 3;
1224  r = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
1225
1226  //Copy
1227  for (i=0; i<3; i++) {
1228    ((double *) (r -> data))[i] = ((double *) (q -> data))[i];
1229  }
1230
1231  //Get normal and direction
1232  n1 = ((double *) normal -> data)[0];
1233  n2 = ((double *) normal -> data)[1];
1234  if (direction == -1) n2 = -n2;
1235
1236  //Rotate
1237  _rotate((double *) r -> data, n1, n2);
1238
1239  //Release numeric arrays
1240  Py_DECREF(q);
1241  Py_DECREF(normal);
1242
1243  //return result using PyArray to avoid memory leak
1244  return PyArray_Return(r);
1245}
1246
1247PyObject *compute_fluxes_ext_central(PyObject *self, PyObject *args) {
1248  /*Compute all fluxes and the timestep suitable for all volumes
1249    in domain.
1250
1251    Compute total flux for each conserved quantity using "flux_function_central"
1252
1253    Fluxes across each edge are scaled by edgelengths and summed up
1254    Resulting flux is then scaled by area and stored in
1255    explicit_update for each of the three conserved quantities
1256    stage, xmomentum and ymomentum
1257
1258    The maximal allowable speed computed by the flux_function for each volume
1259    is converted to a timestep that must not be exceeded. The minimum of
1260    those is computed as the next overall timestep.
1261
1262    Python call:
1263    domain.timestep = compute_fluxes(timestep,
1264                                     domain.epsilon,
1265                                     domain.H0,
1266                                     domain.g,
1267                                     domain.neighbours,
1268                                     domain.neighbour_edges,
1269                                     domain.normals,
1270                                     domain.edgelengths,
1271                                     domain.radii,
1272                                     domain.areas,
1273                                     tri_full_flag,
1274                                     Stage.edge_values,
1275                                     Xmom.edge_values,
1276                                     Ymom.edge_values,
1277                                     Bed.edge_values,
1278                                     Stage.boundary_values,
1279                                     Xmom.boundary_values,
1280                                     Ymom.boundary_values,
1281                                     Stage.explicit_update,
1282                                     Xmom.explicit_update,
1283                                     Ymom.explicit_update,
1284                                     already_computed_flux,
1285                                     optimise_dry_cells)                                     
1286
1287
1288    Post conditions:
1289      domain.explicit_update is reset to computed flux values
1290      domain.timestep is set to the largest step satisfying all volumes.
1291
1292
1293  */
1294
1295
1296  PyArrayObject *neighbours, *neighbour_edges,
1297    *normals, *edgelengths, *radii, *areas,
1298    *tri_full_flag,
1299    *stage_edge_values,
1300    *xmom_edge_values,
1301    *ymom_edge_values,
1302    *bed_edge_values,
1303    *stage_boundary_values,
1304    *xmom_boundary_values,
1305    *ymom_boundary_values,
1306    *stage_explicit_update,
1307    *xmom_explicit_update,
1308    *ymom_explicit_update,
1309    *already_computed_flux, //Tracks whether the flux across an edge has already been computed
1310    *max_speed_array; //Keeps track of max speeds for each triangle
1311
1312
1313  // Local variables
1314  double timestep, max_speed, epsilon, g, H0, length, area;
1315  int optimise_dry_cells=0; // Optimisation flag 
1316  double normal[2], ql[3], qr[3], zl, zr;
1317  double edgeflux[3]; // Work array for summing up fluxes
1318
1319  int number_of_elements, k, i, m, n; //, j, computation_needed;
1320
1321  int ki, nm=0, ki2; // Index shorthands
1322  static long call=1; // Static local variable flagging already computed flux
1323
1324
1325  // Convert Python arguments to C
1326  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOi",
1327                        &timestep,
1328                        &epsilon,
1329                        &H0,
1330                        &g,
1331                        &neighbours,
1332                        &neighbour_edges,
1333                        &normals,
1334                        &edgelengths, &radii, &areas,
1335                        &tri_full_flag,
1336                        &stage_edge_values,
1337                        &xmom_edge_values,
1338                        &ymom_edge_values,
1339                        &bed_edge_values,
1340                        &stage_boundary_values,
1341                        &xmom_boundary_values,
1342                        &ymom_boundary_values,
1343                        &stage_explicit_update,
1344                        &xmom_explicit_update,
1345                        &ymom_explicit_update,
1346                        &already_computed_flux,
1347                        &max_speed_array,
1348                        &optimise_dry_cells)) {
1349    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
1350    return NULL;
1351  }
1352 
1353 
1354  number_of_elements = stage_edge_values -> dimensions[0];
1355
1356  call++; // Flag 'id' of flux calculation for this timestep
1357 
1358  // Set explicit_update to zero for all conserved_quantities.
1359  // This assumes compute_fluxes called before forcing terms
1360  for (k=0; k<number_of_elements; k++) {
1361    ((double *) stage_explicit_update -> data)[k]=0.0;
1362    ((double *) xmom_explicit_update -> data)[k]=0.0;
1363    ((double *) ymom_explicit_update -> data)[k]=0.0; 
1364  }
1365 
1366  // For all triangles
1367  for (k=0; k<number_of_elements; k++) {
1368 
1369    // Loop through neighbours and compute edge flux for each 
1370    for (i=0; i<3; i++) {
1371      ki = k*3+i; // Linear index (triangle k, edge i)
1372     
1373      if (((long *) already_computed_flux->data)[ki] == call)
1374        // We've already computed the flux across this edge
1375        continue;
1376       
1377       
1378      ql[0] = ((double *) stage_edge_values -> data)[ki];
1379      ql[1] = ((double *) xmom_edge_values -> data)[ki];
1380      ql[2] = ((double *) ymom_edge_values -> data)[ki];
1381      zl =    ((double *) bed_edge_values -> data)[ki];
1382
1383      // Quantities at neighbour on nearest face
1384      n = ((long *) neighbours -> data)[ki];
1385      if (n < 0) {
1386        m = -n-1; // Convert negative flag to index
1387       
1388        qr[0] = ((double *) stage_boundary_values -> data)[m];
1389        qr[1] = ((double *) xmom_boundary_values -> data)[m];
1390        qr[2] = ((double *) ymom_boundary_values -> data)[m];
1391        zr = zl; //Extend bed elevation to boundary
1392      } else {
1393        m = ((long *) neighbour_edges -> data)[ki];
1394        nm = n*3+m; // Linear index (triangle n, edge m)
1395       
1396        qr[0] = ((double *) stage_edge_values -> data)[nm];
1397        qr[1] = ((double *) xmom_edge_values -> data)[nm];
1398        qr[2] = ((double *) ymom_edge_values -> data)[nm];
1399        zr =    ((double *) bed_edge_values -> data)[nm];
1400      }
1401     
1402     
1403      if (optimise_dry_cells) {     
1404        // Check if flux calculation is necessary across this edge
1405        // This check will exclude dry cells.
1406        // This will also optimise cases where zl != zr as
1407        // long as both are dry
1408
1409        if ( fabs(ql[0] - zl) < epsilon && 
1410             fabs(qr[0] - zr) < epsilon ) {
1411          // Cell boundary is dry
1412         
1413          ((long *) already_computed_flux -> data)[ki] = call; // #k Done       
1414          if (n>=0)
1415            ((long *) already_computed_flux -> data)[nm] = call; // #n Done
1416       
1417          max_speed = 0.0;
1418          continue;
1419        }
1420      }
1421     
1422
1423     
1424           
1425      // Outward pointing normal vector (domain.normals[k, 2*i:2*i+2])
1426      ki2 = 2*ki; //k*6 + i*2
1427      normal[0] = ((double *) normals -> data)[ki2];
1428      normal[1] = ((double *) normals -> data)[ki2+1];
1429     
1430      // Edge flux computation (triangle k, edge i)
1431      flux_function_central(ql, qr, zl, zr,
1432                            normal[0], normal[1],
1433                            epsilon, H0, g,
1434                            edgeflux, &max_speed);
1435     
1436     
1437      // Multiply edgeflux by edgelength
1438      length = ((double *) edgelengths -> data)[ki];
1439      edgeflux[0] *= length;           
1440      edgeflux[1] *= length;           
1441      edgeflux[2] *= length;                       
1442     
1443     
1444      // Update triangle k with flux from edge i
1445      ((double *) stage_explicit_update -> data)[k] -= edgeflux[0];
1446      ((double *) xmom_explicit_update -> data)[k] -= edgeflux[1];
1447      ((double *) ymom_explicit_update -> data)[k] -= edgeflux[2];
1448     
1449      ((long *) already_computed_flux -> data)[ki] = call; // #k Done
1450     
1451     
1452      // Update neighbour n with same flux but reversed sign
1453      if (n>=0){
1454        ((double *) stage_explicit_update -> data)[n] += edgeflux[0];
1455        ((double *) xmom_explicit_update -> data)[n] += edgeflux[1];
1456        ((double *) ymom_explicit_update -> data)[n] += edgeflux[2];
1457       
1458        ((long *) already_computed_flux -> data)[nm] = call; // #n Done
1459      }
1460     
1461     
1462      // Update timestep based on edge i and possibly neighbour n
1463      if ( ((long *) tri_full_flag -> data)[k] == 1) {
1464        if (max_speed > epsilon) {
1465          timestep = min(timestep, ((double *) radii -> data)[k]/max_speed);
1466          if (n>=0) 
1467            timestep = min(timestep, ((double *) radii -> data)[n]/max_speed);
1468        }
1469      }
1470     
1471    } // End edge i
1472   
1473   
1474    // Normalise triangle k by area and store for when all conserved
1475    // quantities get updated
1476    area = ((double *) areas -> data)[k];
1477    ((double *) stage_explicit_update -> data)[k] /= area;
1478    ((double *) xmom_explicit_update -> data)[k] /= area;
1479    ((double *) ymom_explicit_update -> data)[k] /= area;
1480   
1481   
1482    // Keep track of maximal speeds
1483    ((double *) max_speed_array -> data)[k] = max_speed;   
1484   
1485  } // End triangle k
1486 
1487  return Py_BuildValue("d", timestep);
1488}
1489
1490
1491PyObject *compute_fluxes_ext_kinetic(PyObject *self, PyObject *args) {
1492  /*Compute all fluxes and the timestep suitable for all volumes
1493    in domain.
1494
1495    Compute total flux for each conserved quantity using "flux_function_central"
1496
1497    Fluxes across each edge are scaled by edgelengths and summed up
1498    Resulting flux is then scaled by area and stored in
1499    explicit_update for each of the three conserved quantities
1500    stage, xmomentum and ymomentum
1501
1502    The maximal allowable speed computed by the flux_function for each volume
1503    is converted to a timestep that must not be exceeded. The minimum of
1504    those is computed as the next overall timestep.
1505
1506    Python call:
1507    domain.timestep = compute_fluxes(timestep,
1508                                     domain.epsilon,
1509                                     domain.H0,
1510                                     domain.g,
1511                                     domain.neighbours,
1512                                     domain.neighbour_edges,
1513                                     domain.normals,
1514                                     domain.edgelengths,
1515                                     domain.radii,
1516                                     domain.areas,
1517                                     Stage.edge_values,
1518                                     Xmom.edge_values,
1519                                     Ymom.edge_values,
1520                                     Bed.edge_values,
1521                                     Stage.boundary_values,
1522                                     Xmom.boundary_values,
1523                                     Ymom.boundary_values,
1524                                     Stage.explicit_update,
1525                                     Xmom.explicit_update,
1526                                     Ymom.explicit_update,
1527                                     already_computed_flux)
1528
1529
1530    Post conditions:
1531      domain.explicit_update is reset to computed flux values
1532      domain.timestep is set to the largest step satisfying all volumes.
1533
1534
1535  */
1536
1537
1538  PyArrayObject *neighbours, *neighbour_edges,
1539    *normals, *edgelengths, *radii, *areas,
1540    *tri_full_flag,
1541    *stage_edge_values,
1542    *xmom_edge_values,
1543    *ymom_edge_values,
1544    *bed_edge_values,
1545    *stage_boundary_values,
1546    *xmom_boundary_values,
1547    *ymom_boundary_values,
1548    *stage_explicit_update,
1549    *xmom_explicit_update,
1550    *ymom_explicit_update,
1551    *already_computed_flux;//tracks whether the flux across an edge has already been computed
1552
1553
1554  //Local variables
1555  double timestep, max_speed, epsilon, g, H0;
1556  double normal[2], ql[3], qr[3], zl, zr;
1557  double edgeflux[3]; //Work arrays for summing up fluxes
1558
1559  int number_of_elements, k, i, m, n;
1560  int ki, nm=0, ki2; //Index shorthands
1561  static long call=1;
1562
1563
1564  // Convert Python arguments to C
1565  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOO",
1566                        &timestep,
1567                        &epsilon,
1568                        &H0,
1569                        &g,
1570                        &neighbours,
1571                        &neighbour_edges,
1572                        &normals,
1573                        &edgelengths, &radii, &areas,
1574                        &tri_full_flag,
1575                        &stage_edge_values,
1576                        &xmom_edge_values,
1577                        &ymom_edge_values,
1578                        &bed_edge_values,
1579                        &stage_boundary_values,
1580                        &xmom_boundary_values,
1581                        &ymom_boundary_values,
1582                        &stage_explicit_update,
1583                        &xmom_explicit_update,
1584                        &ymom_explicit_update,
1585                        &already_computed_flux)) {
1586    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
1587    return NULL;
1588  }
1589  number_of_elements = stage_edge_values -> dimensions[0];
1590  call++;//a static local variable to which already_computed_flux is compared
1591  //set explicit_update to zero for all conserved_quantities.
1592  //This assumes compute_fluxes called before forcing terms
1593  for (k=0; k<number_of_elements; k++) {
1594    ((double *) stage_explicit_update -> data)[k]=0.0;
1595    ((double *) xmom_explicit_update -> data)[k]=0.0;
1596    ((double *) ymom_explicit_update -> data)[k]=0.0;
1597  }
1598  //Loop through neighbours and compute edge flux for each
1599  for (k=0; k<number_of_elements; k++) {
1600    for (i=0; i<3; i++) {
1601      ki = k*3+i;
1602      if (((long *) already_computed_flux->data)[ki]==call)//we've already computed the flux across this edge
1603        continue;
1604      ql[0] = ((double *) stage_edge_values -> data)[ki];
1605      ql[1] = ((double *) xmom_edge_values -> data)[ki];
1606      ql[2] = ((double *) ymom_edge_values -> data)[ki];
1607      zl =    ((double *) bed_edge_values -> data)[ki];
1608
1609      //Quantities at neighbour on nearest face
1610      n = ((long *) neighbours -> data)[ki];
1611      if (n < 0) {
1612        m = -n-1; //Convert negative flag to index
1613        qr[0] = ((double *) stage_boundary_values -> data)[m];
1614        qr[1] = ((double *) xmom_boundary_values -> data)[m];
1615        qr[2] = ((double *) ymom_boundary_values -> data)[m];
1616        zr = zl; //Extend bed elevation to boundary
1617      } else {
1618        m = ((long *) neighbour_edges -> data)[ki];
1619        nm = n*3+m;
1620        qr[0] = ((double *) stage_edge_values -> data)[nm];
1621        qr[1] = ((double *) xmom_edge_values -> data)[nm];
1622        qr[2] = ((double *) ymom_edge_values -> data)[nm];
1623        zr =    ((double *) bed_edge_values -> data)[nm];
1624      }
1625      // Outward pointing normal vector
1626      // normal = domain.normals[k, 2*i:2*i+2]
1627      ki2 = 2*ki; //k*6 + i*2
1628      normal[0] = ((double *) normals -> data)[ki2];
1629      normal[1] = ((double *) normals -> data)[ki2+1];
1630      //Edge flux computation
1631      flux_function_kinetic(ql, qr, zl, zr,
1632                    normal[0], normal[1],
1633                    epsilon, H0, g,
1634                    edgeflux, &max_speed);
1635      //update triangle k
1636      ((long *) already_computed_flux->data)[ki]=call;
1637      ((double *) stage_explicit_update -> data)[k] -= edgeflux[0]*((double *) edgelengths -> data)[ki];
1638      ((double *) xmom_explicit_update -> data)[k] -= edgeflux[1]*((double *) edgelengths -> data)[ki];
1639      ((double *) ymom_explicit_update -> data)[k] -= edgeflux[2]*((double *) edgelengths -> data)[ki];
1640      //update the neighbour n
1641      if (n>=0){
1642        ((long *) already_computed_flux->data)[nm]=call;
1643        ((double *) stage_explicit_update -> data)[n] += edgeflux[0]*((double *) edgelengths -> data)[nm];
1644        ((double *) xmom_explicit_update -> data)[n] += edgeflux[1]*((double *) edgelengths -> data)[nm];
1645        ((double *) ymom_explicit_update -> data)[n] += edgeflux[2]*((double *) edgelengths -> data)[nm];
1646      }
1647      ///for (j=0; j<3; j++) {
1648        ///flux[j] -= edgeflux[j]*((double *) edgelengths -> data)[ki];
1649        ///}
1650        //Update timestep
1651        //timestep = min(timestep, domain.radii[k]/max_speed)
1652        //FIXME: SR Add parameter for CFL condition
1653    if ( ((long *) tri_full_flag -> data)[k] == 1) {
1654            if (max_speed > epsilon) {
1655                timestep = min(timestep, ((double *) radii -> data)[k]/max_speed);
1656                //maxspeed in flux_function is calculated as max(|u+a|,|u-a|)
1657                if (n>=0)
1658                    timestep = min(timestep, ((double *) radii -> data)[n]/max_speed);
1659            }
1660    }
1661    } // end for i
1662    //Normalise by area and store for when all conserved
1663    //quantities get updated
1664    ((double *) stage_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
1665    ((double *) xmom_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
1666    ((double *) ymom_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
1667  } //end for k
1668  return Py_BuildValue("d", timestep);
1669}
1670
1671PyObject *protect(PyObject *self, PyObject *args) {
1672  //
1673  //    protect(minimum_allowed_height, maximum_allowed_speed, wc, zc, xmomc, ymomc)
1674
1675
1676  PyArrayObject
1677  *wc,            //Stage at centroids
1678  *zc,            //Elevation at centroids
1679  *xmomc,         //Momentums at centroids
1680  *ymomc;
1681
1682
1683  int N;
1684  double minimum_allowed_height, maximum_allowed_speed, epsilon;
1685
1686  // Convert Python arguments to C
1687  if (!PyArg_ParseTuple(args, "dddOOOO",
1688                        &minimum_allowed_height,
1689                        &maximum_allowed_speed,
1690                        &epsilon,
1691                        &wc, &zc, &xmomc, &ymomc)) {
1692    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: protect could not parse input arguments");
1693    return NULL;
1694  } 
1695
1696  N = wc -> dimensions[0];
1697
1698  _protect(N,
1699           minimum_allowed_height,
1700           maximum_allowed_speed,
1701           epsilon,
1702           (double*) wc -> data,
1703           (double*) zc -> data,
1704           (double*) xmomc -> data,
1705           (double*) ymomc -> data);
1706
1707  return Py_BuildValue("");
1708}
1709
1710
1711
1712PyObject *balance_deep_and_shallow(PyObject *self, PyObject *args) {
1713  //
1714  //    balance_deep_and_shallow(wc, zc, hc, wv, zv, hv,
1715  //                             xmomc, ymomc, xmomv, ymomv)
1716
1717
1718  PyArrayObject
1719    *wc,            //Stage at centroids
1720    *zc,            //Elevation at centroids
1721    *hc,            //Height at centroids
1722    *wv,            //Stage at vertices
1723    *zv,            //Elevation at vertices
1724    *hv,            //Depths at vertices
1725    *hvbar,         //h-Limited depths at vertices
1726    *xmomc,         //Momentums at centroids and vertices
1727    *ymomc,
1728    *xmomv,
1729    *ymomv;
1730 
1731  PyObject *domain, *Tmp;
1732   
1733  double alpha_balance = 2.0;
1734  double H0;
1735
1736  int N, tight_slope_limiters; //, err;
1737
1738  // Convert Python arguments to C
1739  if (!PyArg_ParseTuple(args, "OOOOOOOOOOOO",
1740                        &domain,
1741                        &wc, &zc, &hc,
1742                        &wv, &zv, &hv, &hvbar,
1743                        &xmomc, &ymomc, &xmomv, &ymomv)) {
1744    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not parse input arguments");
1745    return NULL;
1746  } 
1747         
1748  // Pull out parameters
1749  Tmp = PyObject_GetAttrString(domain, "alpha_balance");
1750  if (!Tmp) {
1751    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not obtain object alpha_balance from domain");
1752    return NULL;
1753  } 
1754  alpha_balance = PyFloat_AsDouble(Tmp);
1755  Py_DECREF(Tmp);
1756
1757 
1758  Tmp = PyObject_GetAttrString(domain, "H0");
1759  if (!Tmp) {
1760    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not obtain object H0 from domain");
1761    return NULL;
1762  } 
1763  H0 = PyFloat_AsDouble(Tmp);
1764  Py_DECREF(Tmp);
1765
1766 
1767  Tmp = PyObject_GetAttrString(domain, "tight_slope_limiters");
1768  if (!Tmp) {
1769    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not obtain object tight_slope_limiters from domain");
1770    return NULL;
1771  } 
1772  tight_slope_limiters = PyInt_AsLong(Tmp);
1773  Py_DECREF(Tmp);
1774   
1775
1776 
1777 
1778  N = wc -> dimensions[0];
1779
1780  _balance_deep_and_shallow(N,
1781                            (double*) wc -> data,
1782                            (double*) zc -> data,
1783                            (double*) hc -> data,
1784                            (double*) wv -> data,
1785                            (double*) zv -> data,
1786                            (double*) hv -> data,
1787                            (double*) hvbar -> data,
1788                            (double*) xmomc -> data,
1789                            (double*) ymomc -> data,
1790                            (double*) xmomv -> data,
1791                            (double*) ymomv -> data,
1792                            H0,
1793                            (int) tight_slope_limiters,
1794                            alpha_balance);
1795
1796
1797  return Py_BuildValue("");
1798}
1799
1800
1801
1802PyObject *h_limiter(PyObject *self, PyObject *args) {
1803
1804  PyObject *domain, *Tmp;
1805  PyArrayObject
1806    *hv, *hc, //Depth at vertices and centroids
1807    *hvbar,   //Limited depth at vertices (return values)
1808    *neighbours;
1809
1810  int k, i, n, N, k3;
1811  int dimensions[2];
1812  double beta_h; //Safety factor (see config.py)
1813  double *hmin, *hmax, hn;
1814
1815  // Convert Python arguments to C
1816  if (!PyArg_ParseTuple(args, "OOO", &domain, &hc, &hv)) {
1817    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: h_limiter could not parse input arguments");
1818    return NULL;
1819  } 
1820 
1821  neighbours = get_consecutive_array(domain, "neighbours");
1822
1823  //Get safety factor beta_h
1824  Tmp = PyObject_GetAttrString(domain, "beta_h");
1825  if (!Tmp) {
1826    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: h_limiter could not obtain object beta_h from domain");
1827    return NULL;
1828  } 
1829  beta_h = PyFloat_AsDouble(Tmp);
1830
1831  Py_DECREF(Tmp);
1832
1833  N = hc -> dimensions[0];
1834
1835  //Create hvbar
1836  dimensions[0] = N;
1837  dimensions[1] = 3;
1838  hvbar = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
1839
1840
1841  //Find min and max of this and neighbour's centroid values
1842  hmin = malloc(N * sizeof(double));
1843  hmax = malloc(N * sizeof(double));
1844  for (k=0; k<N; k++) {
1845    k3=k*3;
1846
1847    hmin[k] = ((double*) hc -> data)[k];
1848    hmax[k] = hmin[k];
1849
1850    for (i=0; i<3; i++) {
1851      n = ((long*) neighbours -> data)[k3+i];
1852
1853      //Initialise hvbar with values from hv
1854      ((double*) hvbar -> data)[k3+i] = ((double*) hv -> data)[k3+i];
1855
1856      if (n >= 0) {
1857        hn = ((double*) hc -> data)[n]; //Neighbour's centroid value
1858
1859        hmin[k] = min(hmin[k], hn);
1860        hmax[k] = max(hmax[k], hn);
1861      }
1862    }
1863  }
1864
1865  // Call underlying standard routine
1866  _limit(N, beta_h, (double*) hc -> data, (double*) hvbar -> data, hmin, hmax);
1867
1868  // // //Py_DECREF(domain); //FIXME: NEcessary?
1869  free(hmin);
1870  free(hmax);
1871
1872  //return result using PyArray to avoid memory leak
1873  return PyArray_Return(hvbar);
1874  //return Py_BuildValue("");
1875}
1876
1877PyObject *h_limiter_sw(PyObject *self, PyObject *args) {
1878  //a faster version of h_limiter above
1879  PyObject *domain, *Tmp;
1880  PyArrayObject
1881    *hv, *hc, //Depth at vertices and centroids
1882    *hvbar,   //Limited depth at vertices (return values)
1883    *neighbours;
1884
1885  int k, i, N, k3,k0,k1,k2;
1886  int dimensions[2];
1887  double beta_h; //Safety factor (see config.py)
1888  double hmin, hmax, dh[3];
1889// Convert Python arguments to C
1890  if (!PyArg_ParseTuple(args, "OOO", &domain, &hc, &hv)) {
1891    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: h_limiter_sw could not parse input arguments");
1892    return NULL;
1893  } 
1894  neighbours = get_consecutive_array(domain, "neighbours");
1895
1896  //Get safety factor beta_h
1897  Tmp = PyObject_GetAttrString(domain, "beta_h");
1898  if (!Tmp) {
1899    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: h_limiter_sw could not obtain object beta_h from domain");
1900    return NULL;
1901  } 
1902  beta_h = PyFloat_AsDouble(Tmp);
1903
1904  Py_DECREF(Tmp);
1905
1906  N = hc -> dimensions[0];
1907
1908  //Create hvbar
1909  dimensions[0] = N;
1910  dimensions[1] = 3;
1911  hvbar = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
1912  for (k=0;k<N;k++){
1913    k3=k*3;
1914    //get the ids of the neighbours
1915    k0 = ((long*) neighbours -> data)[k3];
1916    k1 = ((long*) neighbours -> data)[k3+1];
1917    k2 = ((long*) neighbours -> data)[k3+2];
1918    //set hvbar provisionally
1919    for (i=0;i<3;i++){
1920      ((double*) hvbar -> data)[k3+i] = ((double*) hv -> data)[k3+i];
1921      dh[i]=((double*) hvbar -> data)[k3+i]-((double*) hc -> data)[k];
1922    }
1923    hmin=((double*) hc -> data)[k];
1924    hmax=hmin;
1925    if (k0>=0){
1926      hmin=min(hmin,((double*) hc -> data)[k0]);
1927      hmax=max(hmax,((double*) hc -> data)[k0]);
1928    }
1929    if (k1>=0){
1930      hmin=min(hmin,((double*) hc -> data)[k1]);
1931      hmax=max(hmax,((double*) hc -> data)[k1]);
1932    }
1933    if (k2>=0){
1934      hmin=min(hmin,((double*) hc -> data)[k2]);
1935      hmax=max(hmax,((double*) hc -> data)[k2]);
1936    }
1937    hmin-=((double*) hc -> data)[k];
1938    hmax-=((double*) hc -> data)[k];
1939    limit_gradient(dh,hmin,hmax,beta_h);
1940    for (i=0;i<3;i++)
1941      ((double*) hvbar -> data)[k3+i] = ((double*) hc -> data)[k]+dh[i];
1942  }
1943  return PyArray_Return(hvbar);
1944}
1945
1946PyObject *assign_windfield_values(PyObject *self, PyObject *args) {
1947  //
1948  //      assign_windfield_values(xmom_update, ymom_update,
1949  //                              s_vec, phi_vec, self.const)
1950
1951
1952
1953  PyArrayObject   //(one element per triangle)
1954  *s_vec,         //Speeds
1955  *phi_vec,       //Bearings
1956  *xmom_update,   //Momentum updates
1957  *ymom_update;
1958
1959
1960  int N;
1961  double cw;
1962
1963  // Convert Python arguments to C
1964  if (!PyArg_ParseTuple(args, "OOOOd",
1965                        &xmom_update,
1966                        &ymom_update,
1967                        &s_vec, &phi_vec,
1968                        &cw)) {
1969    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: assign_windfield_values could not parse input arguments");
1970    return NULL;
1971  } 
1972                       
1973
1974  N = xmom_update -> dimensions[0];
1975
1976  _assign_wind_field_values(N,
1977           (double*) xmom_update -> data,
1978           (double*) ymom_update -> data,
1979           (double*) s_vec -> data,
1980           (double*) phi_vec -> data,
1981           cw);
1982
1983  return Py_BuildValue("");
1984}
1985
1986
1987
1988
1989//////////////////////////////////////////
1990// Method table for python module
1991static struct PyMethodDef MethodTable[] = {
1992  /* The cast of the function is necessary since PyCFunction values
1993   * only take two PyObject* parameters, and rotate() takes
1994   * three.
1995   */
1996
1997  {"rotate", (PyCFunction)rotate, METH_VARARGS | METH_KEYWORDS, "Print out"},
1998  {"extrapolate_second_order_sw", extrapolate_second_order_sw, METH_VARARGS, "Print out"},
1999  {"compute_fluxes_ext_central", compute_fluxes_ext_central, METH_VARARGS, "Print out"},
2000  {"compute_fluxes_ext_kinetic", compute_fluxes_ext_kinetic, METH_VARARGS, "Print out"},
2001  {"gravity", gravity, METH_VARARGS, "Print out"},
2002  {"manning_friction", manning_friction, METH_VARARGS, "Print out"},
2003  {"balance_deep_and_shallow", balance_deep_and_shallow,
2004   METH_VARARGS, "Print out"},
2005  {"h_limiter", h_limiter,
2006   METH_VARARGS, "Print out"},
2007  {"h_limiter_sw", h_limiter_sw,
2008   METH_VARARGS, "Print out"},
2009  {"protect", protect, METH_VARARGS | METH_KEYWORDS, "Print out"},
2010  {"assign_windfield_values", assign_windfield_values,
2011   METH_VARARGS | METH_KEYWORDS, "Print out"},
2012  //{"distribute_to_vertices_and_edges",
2013  // distribute_to_vertices_and_edges, METH_VARARGS},
2014  //{"update_conserved_quantities",
2015  // update_conserved_quantities, METH_VARARGS},
2016  //{"set_initialcondition",
2017  // set_initialcondition, METH_VARARGS},
2018  {NULL, NULL}
2019};
2020
2021// Module initialisation
2022void initshallow_water_ext(void){
2023  Py_InitModule("shallow_water_ext", MethodTable);
2024
2025  import_array();     //Necessary for handling of NumPY structures
2026}
Note: See TracBrowser for help on using the repository browser.