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

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

Ole and I seem to have got edge limiteri working.
Set domain.use_edge_limiter = True
to test the method

File size: 72.6 KB
RevLine 
[4376]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
[4733]20// Shared code snippets
[4376]21#include "util_ext.h"
22
[4690]23
[4376]24const double pi = 3.14159265358979;
25
26// Computational function for rotation
27int _rotate(double *q, double n1, double n2) {
28  /*Rotate the momentum component q (q[1], q[2])
29    from x,y coordinates to coordinates based on normal vector (n1, n2).
30
31    Result is returned in array 3x1 r
32    To rotate in opposite direction, call rotate with (q, n1, -n2)
33
34    Contents of q are changed by this function */
35
36
37  double q1, q2;
38
[4727]39  // Shorthands
40  q1 = q[1];  // uh momentum
41  q2 = q[2];  // vh momentum
[4376]42
[4727]43  // Rotate
[4376]44  q[1] =  n1*q1 + n2*q2;
45  q[2] = -n2*q1 + n1*q2;
46
47  return 0;
48}
49
[4710]50int find_qmin_and_qmax(double dq0, double dq1, double dq2, 
51                       double *qmin, double *qmax){
52  // Considering the centroid of an FV triangle and the vertices of its
53  // auxiliary triangle, find
54  // qmin=min(q)-qc and qmax=max(q)-qc,
55  // where min(q) and max(q) are respectively min and max over the
56  // four values (at the centroid of the FV triangle and the auxiliary
57  // triangle vertices),
58  // and qc is the centroid
59  // dq0=q(vertex0)-q(centroid of FV triangle)
60  // dq1=q(vertex1)-q(vertex0)
61  // dq2=q(vertex2)-q(vertex0)
[4729]62 
[4376]63  if (dq0>=0.0){
64    if (dq1>=dq2){
65      if (dq1>=0.0)
66        *qmax=dq0+dq1;
67      else
68        *qmax=dq0;
69      if ((*qmin=dq0+dq2)<0)
[4710]70        ;// qmin is already set to correct value
[4376]71      else
72        *qmin=0.0;
73    }
[4710]74    else{// dq1<dq2
[4376]75      if (dq2>0)
76        *qmax=dq0+dq2;
77      else
78        *qmax=dq0;
79      if ((*qmin=dq0+dq1)<0)
[4710]80        ;// qmin is the correct value
[4376]81      else
82        *qmin=0.0;
83    }
84  }
85  else{//dq0<0
86    if (dq1<=dq2){
87      if (dq1<0.0)
88        *qmin=dq0+dq1;
89      else
90        *qmin=dq0;
91      if ((*qmax=dq0+dq2)>0.0)
[4710]92        ;// qmax is already set to the correct value
[4376]93      else
94        *qmax=0.0;
95    }
[4710]96    else{// dq1>dq2
[4376]97      if (dq2<0.0)
98        *qmin=dq0+dq2;
99      else
100        *qmin=dq0;
101      if ((*qmax=dq0+dq1)>0.0)
[4710]102        ;// qmax is already set to the correct value
[4376]103      else
104        *qmax=0.0;
105    }
106  }
107  return 0;
108}
109
110int limit_gradient(double *dqv, double qmin, double qmax, double beta_w){
[4710]111  // Given provisional jumps dqv from the FV triangle centroid to its
112  // vertices and jumps qmin (qmax) between the centroid of the FV
113  // triangle and the minimum (maximum) of the values at the centroid of
114  // the FV triangle and the auxiliary triangle vertices,
115  // calculate a multiplicative factor phi by which the provisional
116  // vertex jumps are to be limited
[4729]117 
[4376]118  int i;
119  double r=1000.0, r0=1.0, phi=1.0;
[4710]120  static double TINY = 1.0e-100; // to avoid machine accuracy problems.
121  // FIXME: Perhaps use the epsilon used elsewhere.
122 
123  // Any provisional jump with magnitude < TINY does not contribute to
124  // the limiting process.
[4729]125  for (i=0;i<3;i++){
[4376]126    if (dqv[i]<-TINY)
127      r0=qmin/dqv[i];
[4729]128     
[4376]129    if (dqv[i]>TINY)
130      r0=qmax/dqv[i];
[4729]131     
[4376]132    r=min(r0,r);
133  }
[4710]134 
[4376]135  phi=min(r*beta_w,1.0);
136  for (i=0;i<3;i++)
137    dqv[i]=dqv[i]*phi;
[4729]138   
[4376]139  return 0;
140}
141
142
[4815]143double compute_froude_number(double uh,
144                             double h, 
145                             double g,
146                             double epsilon) {
[4677]147                         
[4815]148  // Compute Froude number; v/sqrt(gh)
149 
[4677]150  double froude_number;
151 
152  //Compute Froude number (stability diagnostics)
[4815]153  if (h > epsilon) { 
154    froude_number = uh/sqrt(g*h)/h;
155  } else {
156    froude_number = 0.0;
157    // FIXME (Ole): What should it be when dry??
[4677]158  }
[4815]159 
160  return froude_number;
[4677]161}
162
163
164
[4376]165// Function to obtain speed from momentum and depth.
166// This is used by flux functions
167// Input parameters uh and h may be modified by this function.
168double _compute_speed(double *uh, 
169                      double *h, 
170                      double epsilon, 
171                      double h0) {
172 
173  double u;
[4677]174
[4376]175 
176  if (*h < epsilon) {
177    *h = 0.0;  //Could have been negative
178    u = 0.0;
179  } else {
[4677]180    u = *uh/(*h + h0/ *h);   
[4376]181  }
[4677]182 
[4376]183
[4677]184  // Adjust momentum to be consistent with speed
185  *uh = u * *h;
186 
[4376]187  return u;
188}
189
[4726]190// Innermost flux function (using stage w=z+h)
[4769]191int _flux_function_central(double *q_left, double *q_right,
192                           double z_left, double z_right,
193                           double n1, double n2,
194                           double epsilon, double H0, double g,
195                           double *edgeflux, double *max_speed) {
[4376]196
197  /*Compute fluxes between volumes for the shallow water wave equation
198    cast in terms of the 'stage', w = h+z using
199    the 'central scheme' as described in
200
201    Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
202    Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
203    Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.
204
205    The implemented formula is given in equation (3.15) on page 714
206  */
207
208  int i;
209
210  double w_left, h_left, uh_left, vh_left, u_left;
211  double w_right, h_right, uh_right, vh_right, u_right;
[4677]212  double v_left, v_right; 
[4376]213  double s_min, s_max, soundspeed_left, soundspeed_right;
214  double denom, z;
215
[4727]216  // Workspace (allocate once, use many)
217  static double q_left_rotated[3], q_right_rotated[3], flux_right[3], flux_left[3];
218
[4726]219  double h0 = H0*H0; // This ensures a good balance when h approaches H0.
220                     // But evidence suggests that h0 can be as little as
221                     // epsilon!
[4376]222 
[4687]223  // Copy conserved quantities to protect from modification
[4376]224  for (i=0; i<3; i++) {
[4687]225    q_left_rotated[i] = q_left[i];
226    q_right_rotated[i] = q_right[i];
[4376]227  }
228
[4687]229  // Align x- and y-momentum with x-axis
230  _rotate(q_left_rotated, n1, n2);
231  _rotate(q_right_rotated, n1, n2);
[4376]232
[4687]233  z = (z_left+z_right)/2; // Average elevation values
[4376]234
[4687]235  // Compute speeds in x-direction
236  w_left = q_left_rotated[0];         
[4376]237  h_left = w_left-z;
[4687]238  uh_left = q_left_rotated[1];
[4681]239  u_left = _compute_speed(&uh_left, &h_left, epsilon, h0);
[4376]240
[4687]241  w_right = q_right_rotated[0];
[4376]242  h_right = w_right-z;
[4687]243  uh_right = q_right_rotated[1];
[4681]244  u_right = _compute_speed(&uh_right, &h_right, epsilon, h0); 
[4376]245
[4687]246  // Momentum in y-direction
247  vh_left  = q_left_rotated[2];
248  vh_right = q_right_rotated[2];
[4376]249
[4681]250  // Limit y-momentum if necessary 
251  v_left = _compute_speed(&vh_left, &h_left, epsilon, h0);
252  v_right = _compute_speed(&vh_right, &h_right, epsilon, h0);
[4376]253
[4687]254  // Maximal and minimal wave speeds
[4376]255  soundspeed_left  = sqrt(g*h_left);
256  soundspeed_right = sqrt(g*h_right);
257
258  s_max = max(u_left+soundspeed_left, u_right+soundspeed_right);
259  if (s_max < 0.0) s_max = 0.0;
260
261  s_min = min(u_left-soundspeed_left, u_right-soundspeed_right);
[4726]262  if (s_min > 0.0) s_min = 0.0;
[4376]263
[4685]264 
[4687]265  // Flux formulas
[4376]266  flux_left[0] = u_left*h_left;
267  flux_left[1] = u_left*uh_left + 0.5*g*h_left*h_left;
268  flux_left[2] = u_left*vh_left;
269
270  flux_right[0] = u_right*h_right;
271  flux_right[1] = u_right*uh_right + 0.5*g*h_right*h_right;
272  flux_right[2] = u_right*vh_right;
273
[4685]274 
[4687]275  // Flux computation
[4376]276  denom = s_max-s_min;
[4731]277  if (denom < epsilon) { // FIXME (Ole): Try using H0 here
[4376]278    for (i=0; i<3; i++) edgeflux[i] = 0.0;
279    *max_speed = 0.0;
280  } else {
281    for (i=0; i<3; i++) {
282      edgeflux[i] = s_max*flux_left[i] - s_min*flux_right[i];
[4687]283      edgeflux[i] += s_max*s_min*(q_right_rotated[i]-q_left_rotated[i]);
[4376]284      edgeflux[i] /= denom;
285    }
286
[4687]287    // Maximal wavespeed
[4376]288    *max_speed = max(fabs(s_max), fabs(s_min));
289
[4687]290    // Rotate back
[4376]291    _rotate(edgeflux, n1, -n2);
292  }
[4681]293 
[4376]294  return 0;
295}
296
297double erfcc(double x){
298    double z,t,result;
299
300    z=fabs(x);
301    t=1.0/(1.0+0.5*z);
302    result=t*exp(-z*z-1.26551223+t*(1.00002368+t*(.37409196+
303         t*(.09678418+t*(-.18628806+t*(.27886807+t*(-1.13520398+
304         t*(1.48851587+t*(-.82215223+t*.17087277)))))))));
305    if (x < 0.0) result = 2.0-result;
306
307    return result;
308    }
309
310
311
312// Computational function for flux computation (using stage w=z+h)
313int flux_function_kinetic(double *q_left, double *q_right,
314                  double z_left, double z_right,
315                  double n1, double n2,
316                  double epsilon, double H0, double g,
317                  double *edgeflux, double *max_speed) {
318
319  /*Compute fluxes between volumes for the shallow water wave equation
320    cast in terms of the 'stage', w = h+z using
321    the 'central scheme' as described in
322
323    Zhang et. al., Advances in Water Resources, 26(6), 2003, 635-647.
324  */
325
326  int i;
327
328  double w_left, h_left, uh_left, vh_left, u_left, F_left;
329  double w_right, h_right, uh_right, vh_right, u_right, F_right;
330  double s_min, s_max, soundspeed_left, soundspeed_right;
331  double z;
[4687]332  double q_left_rotated[3], q_right_rotated[3];
[4376]333
334  double h0 = H0*H0; //This ensures a good balance when h approaches H0.
335
336  //Copy conserved quantities to protect from modification
337  for (i=0; i<3; i++) {
[4687]338    q_left_rotated[i] = q_left[i];
339    q_right_rotated[i] = q_right[i];
[4376]340  }
341
342  //Align x- and y-momentum with x-axis
[4687]343  _rotate(q_left_rotated, n1, n2);
344  _rotate(q_right_rotated, n1, n2);
[4376]345
346  z = (z_left+z_right)/2; //Take average of field values
347
348  //Compute speeds in x-direction
[4687]349  w_left = q_left_rotated[0];         
[4376]350  h_left = w_left-z;
[4687]351  uh_left = q_left_rotated[1];
[4376]352  u_left =_compute_speed(&uh_left, &h_left, epsilon, h0);
353
[4687]354  w_right = q_right_rotated[0];
[4376]355  h_right = w_right-z;
[4687]356  uh_right = q_right_rotated[1];
[4376]357  u_right =_compute_speed(&uh_right, &h_right, epsilon, h0); 
358
359
360  //Momentum in y-direction
[4687]361  vh_left  = q_left_rotated[2];
362  vh_right = q_right_rotated[2];
[4376]363
364
365  //Maximal and minimal wave speeds
366  soundspeed_left  = sqrt(g*h_left);
367  soundspeed_right = sqrt(g*h_right);
368
369  s_max = max(u_left+soundspeed_left, u_right+soundspeed_right);
370  if (s_max < 0.0) s_max = 0.0;
371
372  s_min = min(u_left-soundspeed_left, u_right-soundspeed_right);
373  if (s_min > 0.0) s_min = 0.0;
374
375
376  F_left  = 0.0;
377  F_right = 0.0;
378  if (h_left > 0.0) F_left = u_left/sqrt(g*h_left);
379  if (h_right > 0.0) F_right = u_right/sqrt(g*h_right);
380
381  for (i=0; i<3; i++) edgeflux[i] = 0.0;
382  *max_speed = 0.0;
383
384  edgeflux[0] = h_left*u_left/2.0*erfcc(-F_left) +  \
385          h_left*sqrt(g*h_left)/2.0/sqrt(pi)*exp(-(F_left*F_left)) + \
386          h_right*u_right/2.0*erfcc(F_right) -  \
387          h_right*sqrt(g*h_right)/2.0/sqrt(pi)*exp(-(F_right*F_right));
388
389  edgeflux[1] = (h_left*u_left*u_left + g/2.0*h_left*h_left)/2.0*erfcc(-F_left) + \
390          u_left*h_left*sqrt(g*h_left)/2.0/sqrt(pi)*exp(-(F_left*F_left)) + \
391          (h_right*u_right*u_right + g/2.0*h_right*h_right)/2.0*erfcc(F_right) -  \
392          u_right*h_right*sqrt(g*h_right)/2.0/sqrt(pi)*exp(-(F_right*F_right));
393
394  edgeflux[2] = vh_left*u_left/2.0*erfcc(-F_left) + \
395          vh_left*sqrt(g*h_left)/2.0/sqrt(pi)*exp(-(F_left*F_left)) + \
396          vh_right*u_right/2.0*erfcc(F_right) - \
397          vh_right*sqrt(g*h_right)/2.0/sqrt(pi)*exp(-(F_right*F_right));
398
399  //Maximal wavespeed
400  *max_speed = max(fabs(s_max), fabs(s_min));
401
402  //Rotate back
403  _rotate(edgeflux, n1, -n2);
404
405  return 0;
406}
407
408
409
410
411void _manning_friction(double g, double eps, int N,
412                       double* w, double* z,
413                       double* uh, double* vh,
414                       double* eta, double* xmom, double* ymom) {
415
416  int k;
417  double S, h;
418
419  for (k=0; k<N; k++) {
420    if (eta[k] > eps) {
421      h = w[k]-z[k];
422      if (h >= eps) {
423        S = -g * eta[k]*eta[k] * sqrt((uh[k]*uh[k] + vh[k]*vh[k]));
424        S /= pow(h, 7.0/3);      //Expensive (on Ole's home computer)
425        //S /= exp(7.0/3.0*log(h));      //seems to save about 15% over manning_friction
426        //S /= h*h*(1 + h/3.0 - h*h/9.0); //FIXME: Could use a Taylor expansion
427
428
429        //Update momentum
430        xmom[k] += S*uh[k];
431        ymom[k] += S*vh[k];
432      }
433    }
434  }
435}
436
437
438/*
439void _manning_friction_explicit(double g, double eps, int N,
440                       double* w, double* z,
441                       double* uh, double* vh,
442                       double* eta, double* xmom, double* ymom) {
443
444  int k;
445  double S, h;
446
447  for (k=0; k<N; k++) {
448    if (eta[k] > eps) {
449      h = w[k]-z[k];
450      if (h >= eps) {
451        S = -g * eta[k]*eta[k] * sqrt((uh[k]*uh[k] + vh[k]*vh[k]));
452        S /= pow(h, 7.0/3);      //Expensive (on Ole's home computer)
453        //S /= exp(7.0/3.0*log(h));      //seems to save about 15% over manning_friction
454        //S /= h*h*(1 + h/3.0 - h*h/9.0); //FIXME: Could use a Taylor expansion
455
456
457        //Update momentum
458        xmom[k] += S*uh[k];
459        ymom[k] += S*vh[k];
460      }
461    }
462  }
463}
464*/
465
[4823]466
467
468
469double velocity_balance(double uh_i, 
470                        double uh,
471                        double h_i, 
472                        double h, 
473                        double alpha,
474                        double epsilon) {
475  // Find alpha such that speed at the vertex is within one
476  // order of magnitude of the centroid speed   
477
478  // FIXME(Ole): Work in progress
479 
480  double a, b, estimate;
481  double m=10; // One order of magnitude - allow for velocity deviations at vertices
482
483 
484  printf("alpha = %f, uh_i=%f, uh=%f, h_i=%f, h=%f\n",
485         alpha, uh_i, uh, h_i, h);
486     
487 
488 
489   
490  // Shorthands and determine inequality
491  if (fabs(uh) < epsilon) {
492    a = 1.0e10; // Limit
493  } else {
494    a = fabs(uh_i - uh)/fabs(uh);
495  }
496     
497  if (h < epsilon) {
498    b = 1.0e10; // Limit
499  } else {       
500    b = m*fabs(h_i - h)/h;
501  }
502
503  printf("a %f, b %f\n", a, b); 
504
505  if (a > b) {
506    estimate = (m-1)/(a-b);             
507   
508    printf("Alpha %f, estimate %f\n", 
509           alpha, estimate);   
510           
511    if (alpha < estimate) {
512      printf("Adjusting alpha from %f to %f\n", 
513             alpha, estimate);
514      alpha = estimate;
515    }
516  } else {
517 
518    if (h < h_i) {
519      estimate = (m-1)/(a-b);                 
520   
521      printf("Alpha %f, estimate %f\n", 
522             alpha, estimate);   
523           
524      if (alpha < estimate) {
525        printf("Adjusting alpha from %f to %f\n", 
526               alpha, estimate);
527        alpha = estimate;
528      }   
529    }
530    // Always fulfilled as alpha and m-1 are non negative
531  }
532 
533 
534  return alpha;
535}
536
537
[4376]538int _balance_deep_and_shallow(int N,
[4733]539                              double beta_h,
[4376]540                              double* wc,
541                              double* zc,
542                              double* wv,
543                              double* zv,
[4733]544                              double* hvbar, // Retire this
[4376]545                              double* xmomc,
546                              double* ymomc,
547                              double* xmomv,
548                              double* ymomv,
549                              double H0,
[4631]550                              int tight_slope_limiters,
[4376]551                              double alpha_balance) {
552
[4815]553  int k, k3, i, excessive_froude_number=0;
[4823]554
[4733]555  double dz, hmin, alpha, h_diff, hc_k;
[4815]556  double epsilon = 1.0e-6; // FIXME: Temporary measure
557  double g = 9.81; // FIXME: Temporary measure
558  double hv[3], h; // Depths at vertices
559  double Fx, Fy; // Froude numbers
[4825]560  double uc, vc; // Centroid speeds
[4376]561
[4714]562  // Compute linear combination between w-limited stages and
563  // h-limited stages close to the bed elevation.
[4376]564
565  for (k=0; k<N; k++) {
566    // Compute maximal variation in bed elevation
567    // This quantitiy is
568    //     dz = max_i abs(z_i - z_c)
569    // and it is independent of dimension
570    // In the 1d case zc = (z0+z1)/2
571    // In the 2d case zc = (z0+z1+z2)/3
572
573    k3 = 3*k;
[4733]574    hc_k = wc[k] - zc[k]; // Centroid value at triangle k
[4714]575   
[4376]576    dz = 0.0;
[4714]577    if (tight_slope_limiters == 0) {     
578      // FIXME: Try with this one precomputed
579      for (i=0; i<3; i++) {
580        dz = max(dz, fabs(zv[k3+i]-zc[k]));
[4376]581      }
582    }
583
[4815]584    // Calculate depth at vertices (possibly negative here!)
[4733]585    hv[0] = wv[k3] -   zv[k3];
586    hv[1] = wv[k3+1] - zv[k3+1];
587    hv[2] = wv[k3+2] - zv[k3+2];       
588   
[4714]589    // Calculate minimal depth across all three vertices
[4733]590    hmin = min(hv[0], min(hv[1], hv[2]));
[5175]591
592    //if (hmin < 0.0 ) {
593    //  printf("hmin = %f\n",hmin);
594    //}
[4376]595   
[4714]596
597    // Create alpha in [0,1], where alpha==0 means using the h-limited
598    // stage and alpha==1 means using the w-limited stage as
599    // computed by the gradient limiter (both 1st or 2nd order)
[4631]600    if (tight_slope_limiters == 0) {     
[4714]601      // If hmin > dz/alpha_balance then alpha = 1 and the bed will have no
602      // effect
603      // If hmin < 0 then alpha = 0 reverting to constant height above bed.
604      // The parameter alpha_balance==2 by default
[4376]605
606     
607      if (dz > 0.0) {
608        alpha = max( min( alpha_balance*hmin/dz, 1.0), 0.0 );     
609      } else {
[4714]610        alpha = 1.0;  // Flat bed
[4376]611      }
612      //printf("Using old style limiter\n");
613     
614    } else {
615
[4714]616      // Tight Slope Limiter (2007)
[4376]617   
618      // Make alpha as large as possible but still ensure that
[4823]619      // final depth is positive and that velocities at vertices
620      // are controlled
[4376]621   
622      if (hmin < H0) {
623        alpha = 1.0;
624        for (i=0; i<3; i++) {
[4823]625       
626          // FIXME (Ole): Simplify (remove) when (if) hvbar gets retired
[4733]627          if (beta_h > epsilon) {
[4823]628            hc_k = hvbar[k3+i]; // Depth to be used at vertices
[4733]629          }
[4823]630         
631          h_diff = hc_k - hv[i];         
[4376]632          if (h_diff <= 0) {
633            // Deep water triangle is further away from bed than
634            // shallow water (hbar < h). Any alpha will do
635         
636          } else { 
637            // Denominator is positive which means that we need some of the
638            // h-limited stage.
639           
[4823]640            alpha = min(alpha, (hc_k - H0)/h_diff);         
[4376]641          }
642        }
643
644        // Ensure alpha in [0,1]
645        if (alpha>1.0) alpha=1.0;
646        if (alpha<0.0) alpha=0.0;
647       
648      } else {
[4823]649        // Use w-limited stage exclusively in deeper water.
[4376]650        alpha = 1.0;       
651      }
[4823]652
[5175]653      //printf("alpha=%f, tight_slope_limiters=%d\n", alpha, tight_slope_limiters);
654
[4823]655      /*     
656      // Experimental code for controlling velocities at vertices.
657      // Adjust alpha (down towards first order) such that
658      // velocities at vertices remain within one order of magnitude
659      // of those at the centroid.
660     
661      for (i=0; i<3; i++) {     
662     
663        // FIXME (Ole): Simplify (remove) when (if) hvbar gets retired
664        if (beta_h > epsilon) {
665          hc_k = hvbar[k3+i]; // Depth to be used at vertices
666        }
667       
668       
669        alpha = velocity_balance(xmomv[k3+i], xmomc[k],
670                                 hv[i], hc_k, alpha, epsilon);
671                                 
672        alpha = velocity_balance(ymomv[k3+i], ymomc[k],
673                                 hv[i], hc_k, alpha, epsilon);                           
674      }
675      */
[4376]676    }
[4823]677           
[4376]678    //printf("k=%d, hmin=%.2f, dz=%.2f, alpha=%.2f, alpha_balance=%.2f\n",
679    //     k, hmin, dz, alpha, alpha_balance);
680
681    //printf("dz = %.3f, alpha = %.8f\n", dz, alpha);
682
683    //  Let
684    //
685    //    wvi be the w-limited stage (wvi = zvi + hvi)
686    //    wvi- be the h-limited state (wvi- = zvi + hvi-)
687    //
688    //
689    //  where i=0,1,2 denotes the vertex ids
690    //
691    //  Weighted balance between w-limited and h-limited stage is
692    //
693    //    wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi)
694    //
695    //  It follows that the updated wvi is
696    //    wvi := zvi + (1-alpha)*hvi- + alpha*hvi
697    //
698    //   Momentum is balanced between constant and limited
699
[4815]700
701    if (alpha < 1) {     
[4376]702      for (i=0; i<3; i++) {
[4733]703        // FIXME (Ole): Simplify when (if) hvbar gets retired       
704        if (beta_h > epsilon) {   
705          wv[k3+i] = zv[k3+i] + (1-alpha)*hvbar[k3+i] + alpha*hv[i];
706        } else {
707          wv[k3+i] = zv[k3+i] + (1-alpha)*hc_k + alpha*hv[i];   
708        }
[4376]709
[4825]710        // Update momentum at vertices
[5175]711
712
713        // Update momentum at vertices
714        // Update momentum as a linear combination of
715        // xmomc and ymomc (shallow) and momentum
716        // from extrapolator xmomv and ymomv (deep).
717
718
719        //xmomv[k3+i] = (1-alpha)*xmomc[k] + alpha*xmomv[k3+i];
720        //ymomv[k3+i] = (1-alpha)*ymomc[k] + alpha*ymomv[k3+i];
721
722        if (tight_slope_limiters == 1) {
723          // FIXME(Ole): Here's what I think (as of 17 Nov 2007)
[4825]724          // we need to do. Simple and efficient:
725       
726          // Speeds at centroids
727          if (hc_k > epsilon) {
728            uc = xmomc[k]/hc_k;
729            vc = ymomc[k]/hc_k;
730          } else {
731            uc = 0.0;
732            vc = 0.0;
733          }
[4826]734         
[4825]735          // Vertex momenta guaranteed to be consistent with depth guaranteeing
736          // controlled speed
[4826]737          hv[i] = wv[k3+i] - zv[k3+i]; // Recompute (balanced) vertex depth
[4825]738          xmomv[k3+i] = uc*hv[i];
[5175]739          ymomv[k3+i] = vc*hv[i];
[4825]740        } else {
741          // Update momentum as a linear combination of
742          // xmomc and ymomc (shallow) and momentum
743          // from extrapolator xmomv and ymomv (deep).
[5175]744          // FIXME (Ole): Is this really needed? Could we use the above
[4826]745          // instead?
[4825]746         
747          xmomv[k3+i] = (1-alpha)*xmomc[k] + alpha*xmomv[k3+i];
748          ymomv[k3+i] = (1-alpha)*ymomc[k] + alpha*ymomv[k3+i];
749       
750        }
[4376]751      }
752    }
[4815]753       
754
755       
[4826]756    if (0) {  // FIXME(Ole): Disabled while testing balancing of velocities above
757      //if (tight_slope_limiters == 1) {               
[4815]758   
759      // Ensure that the Froude number is kept realistic at vertices
760      // FIXME (Ole): I think it could be used to adjust alpha down
761      // whenever Fr is too large. Possible make sure it doesn't deviate
762      // too much from the value at the centroid. I like this idea!
763     
764      // FIXME (Ole): currently only used with tights_SL
[4825]765     
766      // FIXME (Ole): may not be necessary now
[4815]767
768      excessive_froude_number=0;   
769      for (i=0; i<3; i++) {   
770     
771        // Recalculate depth at vertex i
772        h = wv[k3+i] - zv[k3+i];
773       
774        Fx = compute_froude_number(xmomv[k3+i], h, g, epsilon);
775        Fy = compute_froude_number(ymomv[k3+i], h, g, epsilon);
776       
777        if ( (fabs(Fx) > 100.0) || (fabs(Fy) > 100.0)) {
778          // FIXME: use max_froude - or base it on centroid value of F
779          // printf("Excessive Froude number detected: %f or %f\n", Fx, Fy);
780          excessive_froude_number=1;         
781        }
782      }
783     
784      if (excessive_froude_number) { 
785   
786        // printf("Adjusting momentum to first order.\n");
787        // Go back to first order (alpha = 0) for this triangle
788        for (i=0; i<3; i++) {         
789          xmomv[k3+i] = xmomc[k];
790          ymomv[k3+i] = ymomc[k];
791         
792          if (beta_h > epsilon) {         
793            wv[k3+i] = zv[k3+i] + hvbar[k3+i];
794          } else {
795            wv[k3+i] = zv[k3+i] + hc_k; 
796          }
797        }       
798      }
799    }
[4376]800  }
801  return 0;
802}
803
804
805
806int _protect(int N,
807             double minimum_allowed_height,
808             double maximum_allowed_speed,
809             double epsilon,
810             double* wc,
811             double* zc,
812             double* xmomc,
813             double* ymomc) {
814
815  int k;
816  double hc;
817  double u, v, reduced_speed;
818
819
[4731]820  // Protect against initesimal and negative heights 
821  if (maximum_allowed_speed < epsilon) {
822    for (k=0; k<N; k++) {
823      hc = wc[k] - zc[k];
824
825      if (hc < minimum_allowed_height) {
[4376]826       
[4731]827        // Set momentum to zero and ensure h is non negative
[4376]828        xmomc[k] = 0.0;
829        ymomc[k] = 0.0;
[4731]830        if (hc <= 0.0) wc[k] = zc[k];
831      }
832    }
833  } else {
834   
835    // Protect against initesimal and negative heights
836    for (k=0; k<N; k++) {
837      hc = wc[k] - zc[k];
[4376]838
[4731]839      if (hc < minimum_allowed_height) {
840
841        //New code: Adjust momentum to guarantee speeds are physical
842        //          ensure h is non negative
843             
844        if (hc <= 0.0) {
845                wc[k] = zc[k];
846        xmomc[k] = 0.0;
847        ymomc[k] = 0.0;
848        } else {
849          //Reduce excessive speeds derived from division by small hc
850        //FIXME (Ole): This may be unnecessary with new slope limiters
851        //in effect.
852         
853          u = xmomc[k]/hc;
854          if (fabs(u) > maximum_allowed_speed) {
855            reduced_speed = maximum_allowed_speed * u/fabs(u);
856            //printf("Speed (u) has been reduced from %.3f to %.3f\n",
857            //   u, reduced_speed);
858            xmomc[k] = reduced_speed * hc;
859          }
860
861          v = ymomc[k]/hc;
862          if (fabs(v) > maximum_allowed_speed) {
863            reduced_speed = maximum_allowed_speed * v/fabs(v);
864            //printf("Speed (v) has been reduced from %.3f to %.3f\n",
865            //   v, reduced_speed);
866            ymomc[k] = reduced_speed * hc;
867          }
868        }
[4376]869      }
870    }
871  }
872  return 0;
873}
874
875
876
877int _assign_wind_field_values(int N,
878                              double* xmom_update,
879                              double* ymom_update,
880                              double* s_vec,
881                              double* phi_vec,
882                              double cw) {
883
[4731]884  // Assign windfield values to momentum updates
[4376]885
886  int k;
887  double S, s, phi, u, v;
888
889  for (k=0; k<N; k++) {
890
891    s = s_vec[k];
892    phi = phi_vec[k];
893
894    //Convert to radians
895    phi = phi*pi/180;
896
897    //Compute velocity vector (u, v)
898    u = s*cos(phi);
899    v = s*sin(phi);
900
901    //Compute wind stress
902    S = cw * sqrt(u*u + v*v);
903    xmom_update[k] += S*u;
904    ymom_update[k] += S*v;
905  }
906  return 0;
907}
908
909
910
911///////////////////////////////////////////////////////////////////
912// Gateways to Python
913
[4769]914
915
916PyObject *flux_function_central(PyObject *self, PyObject *args) {
917  //
918  // Gateway to innermost flux function.
919  // This is only used by the unit tests as the c implementation is
920  // normally called by compute_fluxes in this module.
921
922
923  PyArrayObject *normal, *ql, *qr,  *edgeflux;
924  double g, epsilon, max_speed, H0, zl, zr;
925
926  if (!PyArg_ParseTuple(args, "OOOddOddd",
927                        &normal, &ql, &qr, &zl, &zr, &edgeflux,
928                        &epsilon, &g, &H0)) {
929    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: flux_function_central could not parse input arguments");
930    return NULL;
931  }
932
933 
934  _flux_function_central((double*) ql -> data, 
935                         (double*) qr -> data, 
936                         zl, 
937                         zr,                                             
938                         ((double*) normal -> data)[0],                                                                         
939                         ((double*) normal -> data)[1],                 
940                         epsilon, H0, g,
941                         (double*) edgeflux -> data, 
942                         &max_speed);
943 
944  return Py_BuildValue("d", max_speed); 
945}
946
947
948
[4376]949PyObject *gravity(PyObject *self, PyObject *args) {
950  //
951  //  gravity(g, h, v, x, xmom, ymom)
952  //
953
954
955  PyArrayObject *h, *v, *x, *xmom, *ymom;
[4687]956  int k, N, k3, k6;
[4376]957  double g, avg_h, zx, zy;
958  double x0, y0, x1, y1, x2, y2, z0, z1, z2;
[4687]959  //double epsilon;
[4376]960
961  if (!PyArg_ParseTuple(args, "dOOOOO",
962                        &g, &h, &v, &x,
963                        &xmom, &ymom)) {
[4687]964    //&epsilon)) {
[4376]965    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: gravity could not parse input arguments");
966    return NULL;
967  }
968
969  N = h -> dimensions[0];
970  for (k=0; k<N; k++) {
971    k3 = 3*k;  // base index
972
[4687]973    // Get bathymetry
974    z0 = ((double*) v -> data)[k3 + 0];
975    z1 = ((double*) v -> data)[k3 + 1];
976    z2 = ((double*) v -> data)[k3 + 2];
[4376]977
[4687]978    // Optimise for flat bed
979    // Note (Ole): This didn't produce measurable speed up.
980    // Revisit later
981    //if (fabs(z0-z1)<epsilon && fabs(z1-z2)<epsilon) {
982    //  continue;
983    //}
[4376]984
[4687]985    // Get average depth from centroid values
986    avg_h = ((double *) h -> data)[k];
987
988    // Compute bed slope
989    k6 = 6*k;  // base index
990 
[4376]991    x0 = ((double*) x -> data)[k6 + 0];
992    y0 = ((double*) x -> data)[k6 + 1];
993    x1 = ((double*) x -> data)[k6 + 2];
994    y1 = ((double*) x -> data)[k6 + 3];
995    x2 = ((double*) x -> data)[k6 + 4];
996    y2 = ((double*) x -> data)[k6 + 5];
997
998
999    _gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2, &zx, &zy);
1000
[4687]1001    // Update momentum
[4376]1002    ((double*) xmom -> data)[k] += -g*zx*avg_h;
1003    ((double*) ymom -> data)[k] += -g*zy*avg_h;
1004  }
1005
1006  return Py_BuildValue("");
1007}
1008
1009
1010PyObject *manning_friction(PyObject *self, PyObject *args) {
1011  //
1012  // manning_friction(g, eps, h, uh, vh, eta, xmom_update, ymom_update)
1013  //
1014
1015
1016  PyArrayObject *w, *z, *uh, *vh, *eta, *xmom, *ymom;
1017  int N;
1018  double g, eps;
1019
1020  if (!PyArg_ParseTuple(args, "ddOOOOOOO",
1021                        &g, &eps, &w, &z, &uh, &vh, &eta,
1022                        &xmom, &ymom)) {
1023    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: manning_friction could not parse input arguments");
1024    return NULL;
1025  }
1026
1027
1028  N = w -> dimensions[0];
1029  _manning_friction(g, eps, N,
1030                    (double*) w -> data,
1031                    (double*) z -> data,
1032                    (double*) uh -> data,
1033                    (double*) vh -> data,
1034                    (double*) eta -> data,
1035                    (double*) xmom -> data,
1036                    (double*) ymom -> data);
1037
1038  return Py_BuildValue("");
1039}
1040
1041
1042/*
1043PyObject *manning_friction_explicit(PyObject *self, PyObject *args) {
1044  //
1045  // manning_friction_explicit(g, eps, h, uh, vh, eta, xmom_update, ymom_update)
1046  //
1047
1048
1049  PyArrayObject *w, *z, *uh, *vh, *eta, *xmom, *ymom;
1050  int N;
1051  double g, eps;
1052
1053  if (!PyArg_ParseTuple(args, "ddOOOOOOO",
1054                        &g, &eps, &w, &z, &uh, &vh, &eta,
1055                        &xmom, &ymom))
1056    return NULL;
1057
1058  N = w -> dimensions[0];
1059  _manning_friction_explicit(g, eps, N,
1060                    (double*) w -> data,
1061                    (double*) z -> data,
1062                    (double*) uh -> data,
1063                    (double*) vh -> data,
1064                    (double*) eta -> data,
1065                    (double*) xmom -> data,
1066                    (double*) ymom -> data);
1067
1068  return Py_BuildValue("");
1069}
1070*/
1071
[4726]1072
1073
[4729]1074// Computational routine
1075int _extrapolate_second_order_sw(int number_of_elements,
1076                                  double epsilon,
1077                                  double minimum_allowed_height,
1078                                  double beta_w,
1079                                  double beta_w_dry,
1080                                  double beta_uh,
1081                                  double beta_uh_dry,
1082                                  double beta_vh,
1083                                  double beta_vh_dry,
1084                                  long* surrogate_neighbours,
1085                                  long* number_of_boundaries,
1086                                  double* centroid_coordinates,
1087                                  double* stage_centroid_values,
1088                                  double* xmom_centroid_values,
1089                                  double* ymom_centroid_values,
1090                                  double* elevation_centroid_values,
1091                                  double* vertex_coordinates,
1092                                  double* stage_vertex_values,
1093                                  double* xmom_vertex_values,
1094                                  double* ymom_vertex_values,
1095                                  double* elevation_vertex_values,
1096                                  int optimise_dry_cells) {
1097                                 
1098                                 
1099
1100  // Local variables
1101  double a, b; // Gradient vector used to calculate vertex values from centroids
1102  int k,k0,k1,k2,k3,k6,coord_index,i;
1103  double x,y,x0,y0,x1,y1,x2,y2,xv0,yv0,xv1,yv1,xv2,yv2; // Vertices of the auxiliary triangle
1104  double dx1,dx2,dy1,dy2,dxv0,dxv1,dxv2,dyv0,dyv1,dyv2,dq0,dq1,dq2,area2;
1105  double dqv[3], qmin, qmax, hmin, hmax;
[5162]1106  double hc, h0, h1, h2, beta_tmp, hfactor;
[4729]1107
1108       
1109  for (k=0; k<number_of_elements; k++) {
1110    k3=k*3;
1111    k6=k*6;
1112
1113   
1114    if (number_of_boundaries[k]==3){
1115      // No neighbours, set gradient on the triangle to zero
1116     
1117      stage_vertex_values[k3]   = stage_centroid_values[k];
1118      stage_vertex_values[k3+1] = stage_centroid_values[k];
1119      stage_vertex_values[k3+2] = stage_centroid_values[k];
1120      xmom_vertex_values[k3]    = xmom_centroid_values[k];
1121      xmom_vertex_values[k3+1]  = xmom_centroid_values[k];
1122      xmom_vertex_values[k3+2]  = xmom_centroid_values[k];
1123      ymom_vertex_values[k3]    = ymom_centroid_values[k];
1124      ymom_vertex_values[k3+1]  = ymom_centroid_values[k];
1125      ymom_vertex_values[k3+2]  = ymom_centroid_values[k];
1126     
1127      continue;
1128    }
1129    else {
1130      // Triangle k has one or more neighbours.
1131      // Get centroid and vertex coordinates of the triangle
1132     
1133      // Get the vertex coordinates
1134      xv0 = vertex_coordinates[k6];   yv0=vertex_coordinates[k6+1];
1135      xv1 = vertex_coordinates[k6+2]; yv1=vertex_coordinates[k6+3];
1136      xv2 = vertex_coordinates[k6+4]; yv2=vertex_coordinates[k6+5];
1137     
1138      // Get the centroid coordinates
1139      coord_index=2*k;
1140      x=centroid_coordinates[coord_index];
1141      y=centroid_coordinates[coord_index+1];
1142     
1143      // Store x- and y- differentials for the vertices of
1144      // triangle k relative to the centroid
1145      dxv0=xv0-x; dxv1=xv1-x; dxv2=xv2-x;
1146      dyv0=yv0-y; dyv1=yv1-y; dyv2=yv2-y;
1147    }
1148
1149
1150           
1151    if (number_of_boundaries[k]<=1){
1152   
1153      //==============================================
1154      // Number of boundaries <= 1
1155      //==============================================   
1156   
1157   
1158      // If no boundaries, auxiliary triangle is formed
1159      // from the centroids of the three neighbours
1160      // If one boundary, auxiliary triangle is formed
1161      // from this centroid and its two neighbours
1162     
1163      k0=surrogate_neighbours[k3];
1164      k1=surrogate_neighbours[k3+1];
1165      k2=surrogate_neighbours[k3+2];
1166     
1167      // Get the auxiliary triangle's vertex coordinates
1168      // (really the centroids of neighbouring triangles)
1169      coord_index=2*k0;
1170      x0=centroid_coordinates[coord_index];
1171      y0=centroid_coordinates[coord_index+1];
1172     
1173      coord_index=2*k1;
1174      x1=centroid_coordinates[coord_index];
1175      y1=centroid_coordinates[coord_index+1];
1176     
1177      coord_index=2*k2;
1178      x2=centroid_coordinates[coord_index];
1179      y2=centroid_coordinates[coord_index+1];
1180     
1181      // Store x- and y- differentials for the vertices
1182      // of the auxiliary triangle
1183      dx1=x1-x0; dx2=x2-x0;
1184      dy1=y1-y0; dy2=y2-y0;
1185     
1186      // Calculate 2*area of the auxiliary triangle
1187      // The triangle is guaranteed to be counter-clockwise     
1188      area2 = dy2*dx1 - dy1*dx2; 
1189     
1190      // If the mesh is 'weird' near the boundary,
1191      // the triangle might be flat or clockwise:
1192      if (area2<=0) {
1193        PyErr_SetString(PyExc_RuntimeError, 
1194                        "shallow_water_ext.c: negative triangle area encountered");
1195        return -1;
1196      } 
1197     
1198      // Calculate heights of neighbouring cells
1199      hc = stage_centroid_values[k]  - elevation_centroid_values[k];
1200      h0 = stage_centroid_values[k0] - elevation_centroid_values[k0];
1201      h1 = stage_centroid_values[k1] - elevation_centroid_values[k1];
1202      h2 = stage_centroid_values[k2] - elevation_centroid_values[k2];
[5162]1203      hmin = min(min(h0,min(h1,h2)),hc);
1204      //hfactor = hc/(hc + 1.0);
1205
1206      hfactor = 0.0;
[5175]1207      if (hmin > 0.001 ) {
1208          hfactor = (hmin-0.001)/(hmin+0.004);
[5162]1209      }
[4729]1210     
1211      if (optimise_dry_cells) {     
1212        // Check if linear reconstruction is necessary for triangle k
1213        // This check will exclude dry cells.
1214
1215        hmax = max(h0,max(h1,h2));     
1216        if (hmax < epsilon) {
1217          continue;
1218        }
1219      }
1220
1221           
1222      //-----------------------------------
1223      // stage
1224      //-----------------------------------     
1225     
1226      // Calculate the difference between vertex 0 of the auxiliary
1227      // triangle and the centroid of triangle k
1228      dq0=stage_centroid_values[k0]-stage_centroid_values[k];
1229     
1230      // Calculate differentials between the vertices
1231      // of the auxiliary triangle (centroids of neighbouring triangles)
1232      dq1=stage_centroid_values[k1]-stage_centroid_values[k0];
1233      dq2=stage_centroid_values[k2]-stage_centroid_values[k0];
1234     
1235      // Calculate the gradient of stage on the auxiliary triangle
1236      a = dy2*dq1 - dy1*dq2;
1237      a /= area2;
1238      b = dx1*dq2 - dx2*dq1;
1239      b /= area2;
1240     
1241      // Calculate provisional jumps in stage from the centroid
1242      // of triangle k to its vertices, to be limited
1243      dqv[0]=a*dxv0+b*dyv0;
1244      dqv[1]=a*dxv1+b*dyv1;
1245      dqv[2]=a*dxv2+b*dyv2;
1246     
1247      // Now we want to find min and max of the centroid and the
1248      // vertices of the auxiliary triangle and compute jumps
1249      // from the centroid to the min and max
1250      find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
1251     
1252      // Playing with dry wet interface
[5162]1253      //hmin = qmin;
1254      //beta_tmp = beta_w_dry;
1255      //if (hmin>minimum_allowed_height)
1256      beta_tmp = beta_w_dry + (beta_w - beta_w_dry) * hfactor;
[4729]1257       
[5162]1258      //printf("min_alled_height = %f\n",minimum_allowed_height);
1259      //printf("hmin = %f\n",hmin);
1260      //printf("beta_w = %f\n",beta_w);
1261      //printf("beta_tmp = %f\n",beta_tmp);
[4729]1262      // Limit the gradient
1263      limit_gradient(dqv,qmin,qmax,beta_tmp); 
1264     
1265      for (i=0;i<3;i++)
1266        stage_vertex_values[k3+i]=stage_centroid_values[k]+dqv[i];
1267     
1268     
1269      //-----------------------------------
1270      // xmomentum
1271      //-----------------------------------           
1272
1273      // Calculate the difference between vertex 0 of the auxiliary
1274      // triangle and the centroid of triangle k     
1275      dq0=xmom_centroid_values[k0]-xmom_centroid_values[k];
1276     
1277      // Calculate differentials between the vertices
1278      // of the auxiliary triangle
1279      dq1=xmom_centroid_values[k1]-xmom_centroid_values[k0];
1280      dq2=xmom_centroid_values[k2]-xmom_centroid_values[k0];
1281     
1282      // Calculate the gradient of xmom on the auxiliary triangle
1283      a = dy2*dq1 - dy1*dq2;
1284      a /= area2;
1285      b = dx1*dq2 - dx2*dq1;
1286      b /= area2;
1287     
1288      // Calculate provisional jumps in stage from the centroid
1289      // of triangle k to its vertices, to be limited     
1290      dqv[0]=a*dxv0+b*dyv0;
1291      dqv[1]=a*dxv1+b*dyv1;
1292      dqv[2]=a*dxv2+b*dyv2;
1293     
1294      // Now we want to find min and max of the centroid and the
1295      // vertices of the auxiliary triangle and compute jumps
1296      // from the centroid to the min and max
1297      find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
[5162]1298      //beta_tmp = beta_uh;
1299      //if (hmin<minimum_allowed_height)
1300      //beta_tmp = beta_uh_dry;
1301      beta_tmp = beta_uh_dry + (beta_uh - beta_uh_dry) * hfactor;
1302
[4729]1303      // Limit the gradient
1304      limit_gradient(dqv,qmin,qmax,beta_tmp);
1305
1306      for (i=0;i<3;i++)
1307        xmom_vertex_values[k3+i]=xmom_centroid_values[k]+dqv[i];
1308     
1309     
1310      //-----------------------------------
1311      // ymomentum
1312      //-----------------------------------                 
1313
1314      // Calculate the difference between vertex 0 of the auxiliary
1315      // triangle and the centroid of triangle k
1316      dq0=ymom_centroid_values[k0]-ymom_centroid_values[k];
1317     
1318      // Calculate differentials between the vertices
1319      // of the auxiliary triangle
1320      dq1=ymom_centroid_values[k1]-ymom_centroid_values[k0];
1321      dq2=ymom_centroid_values[k2]-ymom_centroid_values[k0];
1322     
1323      // Calculate the gradient of xmom on the auxiliary triangle
1324      a = dy2*dq1 - dy1*dq2;
1325      a /= area2;
1326      b = dx1*dq2 - dx2*dq1;
1327      b /= area2;
1328     
1329      // Calculate provisional jumps in stage from the centroid
1330      // of triangle k to its vertices, to be limited
1331      dqv[0]=a*dxv0+b*dyv0;
1332      dqv[1]=a*dxv1+b*dyv1;
1333      dqv[2]=a*dxv2+b*dyv2;
1334     
1335      // Now we want to find min and max of the centroid and the
1336      // vertices of the auxiliary triangle and compute jumps
1337      // from the centroid to the min and max
1338      find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
1339     
[5162]1340      //beta_tmp = beta_vh;
1341      //
1342      //if (hmin<minimum_allowed_height)
1343      //beta_tmp = beta_vh_dry;
1344      beta_tmp = beta_vh_dry + (beta_vh - beta_vh_dry) * hfactor;       
1345
[4729]1346      // Limit the gradient
1347      limit_gradient(dqv,qmin,qmax,beta_tmp);
1348     
1349      for (i=0;i<3;i++)
1350        ymom_vertex_values[k3+i]=ymom_centroid_values[k]+dqv[i];
1351       
1352    } // End number_of_boundaries <=1
1353    else{
1354
1355      //==============================================
1356      // Number of boundaries == 2
1357      //==============================================       
1358       
1359      // One internal neighbour and gradient is in direction of the neighbour's centroid
1360     
1361      // Find the only internal neighbour
1362      for (k2=k3;k2<k3+3;k2++){
1363        // Find internal neighbour of triangle k     
1364        // k2 indexes the edges of triangle k   
1365     
1366        if (surrogate_neighbours[k2]!=k)
1367          break;
1368      }
1369      if ((k2==k3+3)) {
1370        // If we didn't find an internal neighbour
1371        PyErr_SetString(PyExc_RuntimeError, 
1372                        "shallow_water_ext.c: Internal neighbour not found");     
1373        return -1;
1374      }
1375     
1376      k1=surrogate_neighbours[k2];
1377     
1378      // The coordinates of the triangle are already (x,y).
1379      // Get centroid of the neighbour (x1,y1)
1380      coord_index=2*k1;
1381      x1=centroid_coordinates[coord_index];
1382      y1=centroid_coordinates[coord_index+1];
1383     
1384      // Compute x- and y- distances between the centroid of
1385      // triangle k and that of its neighbour
1386      dx1=x1-x; dy1=y1-y;
1387     
1388      // Set area2 as the square of the distance
1389      area2=dx1*dx1+dy1*dy1;
1390     
1391      // Set dx2=(x1-x0)/((x1-x0)^2+(y1-y0)^2)
1392      // and dy2=(y1-y0)/((x1-x0)^2+(y1-y0)^2) which
1393      // respectively correspond to the x- and y- gradients
1394      // of the conserved quantities
1395      dx2=1.0/area2;
1396      dy2=dx2*dy1;
1397      dx2*=dx1;
1398     
1399     
1400     
1401      //-----------------------------------
1402      // stage
1403      //-----------------------------------           
1404
1405      // Compute differentials
1406      dq1=stage_centroid_values[k1]-stage_centroid_values[k];
1407     
1408      // Calculate the gradient between the centroid of triangle k
1409      // and that of its neighbour
1410      a=dq1*dx2;
1411      b=dq1*dy2;
1412     
1413      // Calculate provisional vertex jumps, to be limited
1414      dqv[0]=a*dxv0+b*dyv0;
1415      dqv[1]=a*dxv1+b*dyv1;
1416      dqv[2]=a*dxv2+b*dyv2;
1417     
1418      // Now limit the jumps
1419      if (dq1>=0.0){
1420        qmin=0.0;
1421        qmax=dq1;
1422      }
1423      else{
1424        qmin=dq1;
1425        qmax=0.0;
1426      }
1427     
1428      // Limit the gradient
1429      limit_gradient(dqv,qmin,qmax,beta_w);
1430     
1431      for (i=0;i<3;i++)
1432        stage_vertex_values[k3+i]=stage_centroid_values[k]+dqv[i];
1433     
1434     
1435      //-----------------------------------
1436      // xmomentum
1437      //-----------------------------------                       
1438     
1439      // Compute differentials
1440      dq1=xmom_centroid_values[k1]-xmom_centroid_values[k];
1441     
1442      // Calculate the gradient between the centroid of triangle k
1443      // and that of its neighbour
1444      a=dq1*dx2;
1445      b=dq1*dy2;
1446     
1447      // Calculate provisional vertex jumps, to be limited
1448      dqv[0]=a*dxv0+b*dyv0;
1449      dqv[1]=a*dxv1+b*dyv1;
1450      dqv[2]=a*dxv2+b*dyv2;
1451     
1452      // Now limit the jumps
1453      if (dq1>=0.0){
1454        qmin=0.0;
1455        qmax=dq1;
1456      }
1457      else{
1458        qmin=dq1;
1459        qmax=0.0;
1460      }
1461     
1462      // Limit the gradient     
1463      limit_gradient(dqv,qmin,qmax,beta_w);
1464     
1465      for (i=0;i<3;i++)
1466        xmom_vertex_values[k3+i]=xmom_centroid_values[k]+dqv[i];
1467     
1468     
1469      //-----------------------------------
1470      // ymomentum
1471      //-----------------------------------                       
1472
1473      // Compute differentials
1474      dq1=ymom_centroid_values[k1]-ymom_centroid_values[k];
1475     
1476      // Calculate the gradient between the centroid of triangle k
1477      // and that of its neighbour
1478      a=dq1*dx2;
1479      b=dq1*dy2;
1480     
1481      // Calculate provisional vertex jumps, to be limited
1482      dqv[0]=a*dxv0+b*dyv0;
1483      dqv[1]=a*dxv1+b*dyv1;
1484      dqv[2]=a*dxv2+b*dyv2;
1485     
1486      // Now limit the jumps
1487      if (dq1>=0.0){
1488        qmin=0.0;
1489        qmax=dq1;
1490      }
1491      else{
1492        qmin=dq1;
1493        qmax=0.0;
1494      }
1495     
1496      // Limit the gradient
1497      limit_gradient(dqv,qmin,qmax,beta_w);
1498     
1499      for (i=0;i<3;i++)
1500        ymom_vertex_values[k3+i]=ymom_centroid_values[k]+dqv[i];
1501       
1502    } // else [number_of_boundaries==2]
1503  } // for k=0 to number_of_elements-1
1504 
1505  return 0;
1506}                       
1507
1508
[4376]1509PyObject *extrapolate_second_order_sw(PyObject *self, PyObject *args) {
[4730]1510  /*Compute the vertex values based on a linear reconstruction
1511    on each triangle
1512   
[4376]1513    These values are calculated as follows:
[4730]1514    1) For each triangle not adjacent to a boundary, we consider the
1515       auxiliary triangle formed by the centroids of its three
1516       neighbours.
1517    2) For each conserved quantity, we integrate around the auxiliary
1518       triangle's boundary the product of the quantity and the outward
1519       normal vector. Dividing by the triangle area gives (a,b), the
1520       average of the vector (q_x,q_y) on the auxiliary triangle.
1521       We suppose that the linear reconstruction on the original
1522       triangle has gradient (a,b).
1523    3) Provisional vertex jumps dqv[0,1,2] are computed and these are
1524       then limited by calling the functions find_qmin_and_qmax and
1525       limit_gradient
[4376]1526
1527    Python call:
1528    extrapolate_second_order_sw(domain.surrogate_neighbours,
1529                                domain.number_of_boundaries
1530                                domain.centroid_coordinates,
1531                                Stage.centroid_values
1532                                Xmom.centroid_values
1533                                Ymom.centroid_values
1534                                domain.vertex_coordinates,
1535                                Stage.vertex_values,
1536                                Xmom.vertex_values,
1537                                Ymom.vertex_values)
1538
1539    Post conditions:
[4730]1540            The vertices of each triangle have values from a
1541            limited linear reconstruction
[4376]1542            based on centroid values
1543
1544  */
1545  PyArrayObject *surrogate_neighbours,
1546    *number_of_boundaries,
1547    *centroid_coordinates,
1548    *stage_centroid_values,
1549    *xmom_centroid_values,
1550    *ymom_centroid_values,
[4978]1551    *elevation_centroid_values,
[4376]1552    *vertex_coordinates,
1553    *stage_vertex_values,
1554    *xmom_vertex_values,
1555    *ymom_vertex_values,
[4978]1556    *elevation_vertex_values;
1557 
1558  PyObject *domain;
[4729]1559
1560 
1561  double beta_w, beta_w_dry, beta_uh, beta_uh_dry, beta_vh, beta_vh_dry;   
1562  double minimum_allowed_height, epsilon;
1563  int optimise_dry_cells, number_of_elements, e;
1564 
1565  // Provisional jumps from centroids to v'tices and safety factor re limiting
1566  // by which these jumps are limited
1567  // Convert Python arguments to C
1568  if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOi",
1569                        &domain,
1570                        &surrogate_neighbours,
1571                        &number_of_boundaries,
1572                        &centroid_coordinates,
1573                        &stage_centroid_values,
1574                        &xmom_centroid_values,
1575                        &ymom_centroid_values,
1576                        &elevation_centroid_values,
1577                        &vertex_coordinates,
1578                        &stage_vertex_values,
1579                        &xmom_vertex_values,
1580                        &ymom_vertex_values,
1581                        &elevation_vertex_values,
1582                        &optimise_dry_cells)) {                 
1583                       
[4730]1584    PyErr_SetString(PyExc_RuntimeError, 
1585                    "Input arguments to extrapolate_second_order_sw failed");
[4729]1586    return NULL;
1587  }
1588
[4730]1589  // Get the safety factor beta_w, set in the config.py file.
1590  // This is used in the limiting process
[4729]1591 
1592
[5162]1593  beta_w                 = get_python_double(domain,"beta_w");
1594  beta_w_dry             = get_python_double(domain,"beta_w_dry");
1595  beta_uh                = get_python_double(domain,"beta_uh");
1596  beta_uh_dry            = get_python_double(domain,"beta_uh_dry");
1597  beta_vh                = get_python_double(domain,"beta_vh");
1598  beta_vh_dry            = get_python_double(domain,"beta_vh_dry"); 
[4729]1599
[5162]1600  minimum_allowed_height = get_python_double(domain,"minimum_allowed_height");
1601  epsilon                = get_python_double(domain,"epsilon");
[4978]1602
1603  number_of_elements = stage_centroid_values -> dimensions[0]; 
1604
[4729]1605  // Call underlying computational routine
1606  e = _extrapolate_second_order_sw(number_of_elements,
1607                                   epsilon,
1608                                   minimum_allowed_height,
1609                                   beta_w,
1610                                   beta_w_dry,
1611                                   beta_uh,
1612                                   beta_uh_dry,
1613                                   beta_vh,
1614                                   beta_vh_dry,
1615                                   (long*) surrogate_neighbours -> data,
1616                                   (long*) number_of_boundaries -> data,
1617                                   (double*) centroid_coordinates -> data,
1618                                   (double*) stage_centroid_values -> data,
1619                                   (double*) xmom_centroid_values -> data,
1620                                   (double*) ymom_centroid_values -> data,
1621                                   (double*) elevation_centroid_values -> data,
1622                                   (double*) vertex_coordinates -> data,
1623                                   (double*) stage_vertex_values -> data,
1624                                   (double*) xmom_vertex_values -> data,
1625                                   (double*) ymom_vertex_values -> data,
1626                                   (double*) elevation_vertex_values -> data,
1627                                   optimise_dry_cells);
1628  if (e == -1) {
1629    // Use error string set inside computational routine
1630    return NULL;
1631  }                               
1632 
1633 
1634  return Py_BuildValue("");
1635 
1636}// extrapolate_second-order_sw
1637
1638
1639
1640
[4376]1641PyObject *rotate(PyObject *self, PyObject *args, PyObject *kwargs) {
1642  //
1643  // r = rotate(q, normal, direction=1)
1644  //
1645  // Where q is assumed to be a Float numeric array of length 3 and
1646  // normal a Float numeric array of length 2.
1647
[4682]1648  // FIXME(Ole): I don't think this is used anymore
[4376]1649
1650  PyObject *Q, *Normal;
1651  PyArrayObject *q, *r, *normal;
1652
1653  static char *argnames[] = {"q", "normal", "direction", NULL};
1654  int dimensions[1], i, direction=1;
1655  double n1, n2;
1656
1657  // Convert Python arguments to C
1658  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OO|i", argnames,
1659                                   &Q, &Normal, &direction)) {
1660    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: rotate could not parse input arguments");
1661    return NULL;
1662  } 
1663
[4769]1664  // Input checks (convert sequences into numeric arrays)
[4376]1665  q = (PyArrayObject *)
1666    PyArray_ContiguousFromObject(Q, PyArray_DOUBLE, 0, 0);
1667  normal = (PyArrayObject *)
1668    PyArray_ContiguousFromObject(Normal, PyArray_DOUBLE, 0, 0);
1669
1670
1671  if (normal -> dimensions[0] != 2) {
1672    PyErr_SetString(PyExc_RuntimeError, "Normal vector must have 2 components");
1673    return NULL;
1674  }
1675
[4769]1676  // Allocate space for return vector r (don't DECREF)
[4376]1677  dimensions[0] = 3;
1678  r = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
1679
[4769]1680  // Copy
[4376]1681  for (i=0; i<3; i++) {
1682    ((double *) (r -> data))[i] = ((double *) (q -> data))[i];
1683  }
1684
[4769]1685  // Get normal and direction
[4376]1686  n1 = ((double *) normal -> data)[0];
1687  n2 = ((double *) normal -> data)[1];
1688  if (direction == -1) n2 = -n2;
1689
[4769]1690  // Rotate
[4376]1691  _rotate((double *) r -> data, n1, n2);
1692
[4769]1693  // Release numeric arrays
[4376]1694  Py_DECREF(q);
1695  Py_DECREF(normal);
1696
[4769]1697  // Return result using PyArray to avoid memory leak
[4376]1698  return PyArray_Return(r);
1699}
1700
[4726]1701
1702// Computational function for flux computation
1703double _compute_fluxes_central(int number_of_elements,
1704                               double timestep,
1705                               double epsilon,
1706                               double H0,
1707                               double g,
1708                               long* neighbours,
1709                               long* neighbour_edges,
1710                               double* normals,
1711                               double* edgelengths, 
1712                               double* radii, 
1713                               double* areas,
1714                               long* tri_full_flag,
1715                               double* stage_edge_values,
1716                               double* xmom_edge_values,
1717                               double* ymom_edge_values,
1718                               double* bed_edge_values,
1719                               double* stage_boundary_values,
1720                               double* xmom_boundary_values,
1721                               double* ymom_boundary_values,
1722                               double* stage_explicit_update,
1723                               double* xmom_explicit_update,
1724                               double* ymom_explicit_update,
1725                               long* already_computed_flux,
1726                               double* max_speed_array,
1727                               int optimise_dry_cells) {
1728                               
1729  // Local variables
[4727]1730  double max_speed, length, area, zl, zr;
1731  int k, i, m, n;
1732  int ki, nm=0, ki2; // Index shorthands
[4726]1733 
[4727]1734  // Workspace (making them static actually made function slightly slower (Ole)) 
1735  double ql[3], qr[3], edgeflux[3]; // Work array for summing up fluxes
[4726]1736
1737  static long call=1; // Static local variable flagging already computed flux
1738                               
1739       
1740  // Start computation
1741  call++; // Flag 'id' of flux calculation for this timestep
1742 
1743  // Set explicit_update to zero for all conserved_quantities.
1744  // This assumes compute_fluxes called before forcing terms
1745  for (k=0; k<number_of_elements; k++) {
1746    stage_explicit_update[k]=0.0;
1747    xmom_explicit_update[k]=0.0;
1748    ymom_explicit_update[k]=0.0; 
1749  }
1750
1751  // For all triangles
1752  for (k=0; k<number_of_elements; k++) {
1753   
1754    // Loop through neighbours and compute edge flux for each 
1755    for (i=0; i<3; i++) {
1756      ki = k*3+i; // Linear index (triangle k, edge i)
1757     
1758      if (already_computed_flux[ki] == call)
1759        // We've already computed the flux across this edge
1760        continue;
1761       
1762       
1763      ql[0] = stage_edge_values[ki];
1764      ql[1] = xmom_edge_values[ki];
1765      ql[2] = ymom_edge_values[ki];
[4727]1766      zl = bed_edge_values[ki];
[4726]1767
1768      // Quantities at neighbour on nearest face
1769      n = neighbours[ki];
1770      if (n < 0) {
1771        m = -n-1; // Convert negative flag to index
1772       
1773        qr[0] = stage_boundary_values[m];
1774        qr[1] = xmom_boundary_values[m];
1775        qr[2] = ymom_boundary_values[m];
1776        zr = zl; // Extend bed elevation to boundary
1777      } else {
1778        m = neighbour_edges[ki];
1779        nm = n*3+m; // Linear index (triangle n, edge m)
1780       
1781        qr[0] = stage_edge_values[nm];
1782        qr[1] = xmom_edge_values[nm];
1783        qr[2] = ymom_edge_values[nm];
1784        zr = bed_edge_values[nm];
1785      }
1786     
1787     
1788      if (optimise_dry_cells) {     
1789        // Check if flux calculation is necessary across this edge
1790        // This check will exclude dry cells.
1791        // This will also optimise cases where zl != zr as
1792        // long as both are dry
1793
1794        if ( fabs(ql[0] - zl) < epsilon && 
1795             fabs(qr[0] - zr) < epsilon ) {
1796          // Cell boundary is dry
1797         
1798          already_computed_flux[ki] = call; // #k Done 
1799          if (n>=0)
1800            already_computed_flux[nm] = call; // #n Done
1801       
1802          max_speed = 0.0;
1803          continue;
1804        }
1805      }
1806   
1807     
1808      // Outward pointing normal vector (domain.normals[k, 2*i:2*i+2])
1809      ki2 = 2*ki; //k*6 + i*2
1810
1811      // Edge flux computation (triangle k, edge i)
[4769]1812      _flux_function_central(ql, qr, zl, zr,
1813                             normals[ki2], normals[ki2+1],
1814                             epsilon, H0, g,
1815                             edgeflux, &max_speed);
[4726]1816     
1817     
1818      // Multiply edgeflux by edgelength
1819      length = edgelengths[ki];
1820      edgeflux[0] *= length;           
1821      edgeflux[1] *= length;           
1822      edgeflux[2] *= length;                       
1823     
1824     
1825      // Update triangle k with flux from edge i
1826      stage_explicit_update[k] -= edgeflux[0];
1827      xmom_explicit_update[k] -= edgeflux[1];
1828      ymom_explicit_update[k] -= edgeflux[2];
1829     
1830      already_computed_flux[ki] = call; // #k Done
1831     
1832     
1833      // Update neighbour n with same flux but reversed sign
1834      if (n>=0) {
1835        stage_explicit_update[n] += edgeflux[0];
1836        xmom_explicit_update[n] += edgeflux[1];
1837        ymom_explicit_update[n] += edgeflux[2];
1838       
1839        already_computed_flux[nm] = call; // #n Done
1840      }
1841
1842
1843      // Update timestep based on edge i and possibly neighbour n
1844      if (tri_full_flag[k] == 1) {
1845        if (max_speed > epsilon) {
1846          timestep = min(timestep, radii[k]/max_speed);
1847          if (n>=0) 
1848            timestep = min(timestep, radii[n]/max_speed);
1849        }
1850      }
1851     
1852    } // End edge i
1853   
1854   
1855    // Normalise triangle k by area and store for when all conserved
1856    // quantities get updated
1857    area = areas[k];
1858    stage_explicit_update[k] /= area;
1859    xmom_explicit_update[k] /= area;
1860    ymom_explicit_update[k] /= area;
1861   
1862   
[4834]1863    // Keep track of maximal speeds
1864    max_speed_array[k] = max_speed;
[4726]1865   
1866  } // End triangle k
1867       
1868                               
1869                               
1870  return timestep;
1871}
1872
1873
[4376]1874PyObject *compute_fluxes_ext_central(PyObject *self, PyObject *args) {
1875  /*Compute all fluxes and the timestep suitable for all volumes
1876    in domain.
1877
1878    Compute total flux for each conserved quantity using "flux_function_central"
1879
1880    Fluxes across each edge are scaled by edgelengths and summed up
1881    Resulting flux is then scaled by area and stored in
1882    explicit_update for each of the three conserved quantities
1883    stage, xmomentum and ymomentum
1884
1885    The maximal allowable speed computed by the flux_function for each volume
1886    is converted to a timestep that must not be exceeded. The minimum of
1887    those is computed as the next overall timestep.
1888
1889    Python call:
1890    domain.timestep = compute_fluxes(timestep,
1891                                     domain.epsilon,
1892                                     domain.H0,
1893                                     domain.g,
1894                                     domain.neighbours,
1895                                     domain.neighbour_edges,
1896                                     domain.normals,
1897                                     domain.edgelengths,
1898                                     domain.radii,
1899                                     domain.areas,
1900                                     tri_full_flag,
1901                                     Stage.edge_values,
1902                                     Xmom.edge_values,
1903                                     Ymom.edge_values,
1904                                     Bed.edge_values,
1905                                     Stage.boundary_values,
1906                                     Xmom.boundary_values,
1907                                     Ymom.boundary_values,
1908                                     Stage.explicit_update,
1909                                     Xmom.explicit_update,
1910                                     Ymom.explicit_update,
[4685]1911                                     already_computed_flux,
1912                                     optimise_dry_cells)                                     
[4376]1913
1914
1915    Post conditions:
1916      domain.explicit_update is reset to computed flux values
1917      domain.timestep is set to the largest step satisfying all volumes.
1918
1919
1920  */
1921
1922
1923  PyArrayObject *neighbours, *neighbour_edges,
1924    *normals, *edgelengths, *radii, *areas,
[4726]1925    *tri_full_flag,
1926    *stage_edge_values,
1927    *xmom_edge_values,
1928    *ymom_edge_values,
1929    *bed_edge_values,
1930    *stage_boundary_values,
1931    *xmom_boundary_values,
1932    *ymom_boundary_values,
1933    *stage_explicit_update,
1934    *xmom_explicit_update,
1935    *ymom_explicit_update,
1936    *already_computed_flux, //Tracks whether the flux across an edge has already been computed
1937    *max_speed_array; //Keeps track of max speeds for each triangle
1938
1939   
1940  double timestep, epsilon, H0, g;
1941  int optimise_dry_cells;
1942   
1943  // Convert Python arguments to C
1944  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOi",
1945                        &timestep,
1946                        &epsilon,
1947                        &H0,
1948                        &g,
1949                        &neighbours,
1950                        &neighbour_edges,
1951                        &normals,
1952                        &edgelengths, &radii, &areas,
1953                        &tri_full_flag,
1954                        &stage_edge_values,
1955                        &xmom_edge_values,
1956                        &ymom_edge_values,
1957                        &bed_edge_values,
1958                        &stage_boundary_values,
1959                        &xmom_boundary_values,
1960                        &ymom_boundary_values,
1961                        &stage_explicit_update,
1962                        &xmom_explicit_update,
1963                        &ymom_explicit_update,
1964                        &already_computed_flux,
1965                        &max_speed_array,
1966                        &optimise_dry_cells)) {
1967    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
1968    return NULL;
1969  }
1970
1971 
1972  int number_of_elements = stage_edge_values -> dimensions[0];
1973
1974  // Call underlying flux computation routine and update
1975  // the explicit update arrays
1976  timestep = _compute_fluxes_central(number_of_elements,
1977                                     timestep,
1978                                     epsilon,
1979                                     H0,
1980                                     g,
[4730]1981                                     (long*) neighbours -> data,
[4726]1982                                     (long*) neighbour_edges -> data,
1983                                     (double*) normals -> data,
1984                                     (double*) edgelengths -> data, 
[4730]1985                                     (double*) radii -> data, 
1986                                     (double*) areas -> data,
[4726]1987                                     (long*) tri_full_flag -> data,
1988                                     (double*) stage_edge_values -> data,
1989                                     (double*) xmom_edge_values -> data,
1990                                     (double*) ymom_edge_values -> data,
1991                                     (double*) bed_edge_values -> data,
1992                                     (double*) stage_boundary_values -> data,
1993                                     (double*) xmom_boundary_values -> data,
1994                                     (double*) ymom_boundary_values -> data,
[4730]1995                                     (double*) stage_explicit_update -> data,
[4726]1996                                     (double*) xmom_explicit_update -> data,
1997                                     (double*) ymom_explicit_update -> data,
1998                                     (long*) already_computed_flux -> data,
1999                                     (double*) max_speed_array -> data,
2000                                     optimise_dry_cells);
2001 
2002  // Return updated flux timestep
2003  return Py_BuildValue("d", timestep);
2004}
2005
2006
2007
[4376]2008PyObject *compute_fluxes_ext_kinetic(PyObject *self, PyObject *args) {
2009  /*Compute all fluxes and the timestep suitable for all volumes
2010    in domain.
2011
2012    Compute total flux for each conserved quantity using "flux_function_central"
2013
2014    Fluxes across each edge are scaled by edgelengths and summed up
2015    Resulting flux is then scaled by area and stored in
2016    explicit_update for each of the three conserved quantities
2017    stage, xmomentum and ymomentum
2018
2019    The maximal allowable speed computed by the flux_function for each volume
2020    is converted to a timestep that must not be exceeded. The minimum of
2021    those is computed as the next overall timestep.
2022
2023    Python call:
2024    domain.timestep = compute_fluxes(timestep,
2025                                     domain.epsilon,
2026                                     domain.H0,
2027                                     domain.g,
2028                                     domain.neighbours,
2029                                     domain.neighbour_edges,
2030                                     domain.normals,
2031                                     domain.edgelengths,
2032                                     domain.radii,
2033                                     domain.areas,
2034                                     Stage.edge_values,
2035                                     Xmom.edge_values,
2036                                     Ymom.edge_values,
2037                                     Bed.edge_values,
2038                                     Stage.boundary_values,
2039                                     Xmom.boundary_values,
2040                                     Ymom.boundary_values,
2041                                     Stage.explicit_update,
2042                                     Xmom.explicit_update,
2043                                     Ymom.explicit_update,
2044                                     already_computed_flux)
2045
2046
2047    Post conditions:
2048      domain.explicit_update is reset to computed flux values
2049      domain.timestep is set to the largest step satisfying all volumes.
2050
2051
2052  */
2053
2054
2055  PyArrayObject *neighbours, *neighbour_edges,
2056    *normals, *edgelengths, *radii, *areas,
2057    *tri_full_flag,
2058    *stage_edge_values,
2059    *xmom_edge_values,
2060    *ymom_edge_values,
2061    *bed_edge_values,
2062    *stage_boundary_values,
2063    *xmom_boundary_values,
2064    *ymom_boundary_values,
2065    *stage_explicit_update,
2066    *xmom_explicit_update,
2067    *ymom_explicit_update,
[4769]2068    *already_computed_flux; // Tracks whether the flux across an edge has already been computed
[4376]2069
2070
[4769]2071  // Local variables
[4376]2072  double timestep, max_speed, epsilon, g, H0;
2073  double normal[2], ql[3], qr[3], zl, zr;
2074  double edgeflux[3]; //Work arrays for summing up fluxes
2075
2076  int number_of_elements, k, i, m, n;
2077  int ki, nm=0, ki2; //Index shorthands
2078  static long call=1;
2079
2080
2081  // Convert Python arguments to C
2082  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOO",
2083                        &timestep,
2084                        &epsilon,
2085                        &H0,
2086                        &g,
2087                        &neighbours,
2088                        &neighbour_edges,
2089                        &normals,
2090                        &edgelengths, &radii, &areas,
2091                        &tri_full_flag,
2092                        &stage_edge_values,
2093                        &xmom_edge_values,
2094                        &ymom_edge_values,
2095                        &bed_edge_values,
2096                        &stage_boundary_values,
2097                        &xmom_boundary_values,
2098                        &ymom_boundary_values,
2099                        &stage_explicit_update,
2100                        &xmom_explicit_update,
2101                        &ymom_explicit_update,
2102                        &already_computed_flux)) {
2103    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
2104    return NULL;
2105  }
2106  number_of_elements = stage_edge_values -> dimensions[0];
2107  call++;//a static local variable to which already_computed_flux is compared
2108  //set explicit_update to zero for all conserved_quantities.
2109  //This assumes compute_fluxes called before forcing terms
2110  for (k=0; k<number_of_elements; k++) {
2111    ((double *) stage_explicit_update -> data)[k]=0.0;
2112    ((double *) xmom_explicit_update -> data)[k]=0.0;
2113    ((double *) ymom_explicit_update -> data)[k]=0.0;
2114  }
2115  //Loop through neighbours and compute edge flux for each
2116  for (k=0; k<number_of_elements; k++) {
2117    for (i=0; i<3; i++) {
2118      ki = k*3+i;
2119      if (((long *) already_computed_flux->data)[ki]==call)//we've already computed the flux across this edge
2120        continue;
2121      ql[0] = ((double *) stage_edge_values -> data)[ki];
2122      ql[1] = ((double *) xmom_edge_values -> data)[ki];
2123      ql[2] = ((double *) ymom_edge_values -> data)[ki];
2124      zl =    ((double *) bed_edge_values -> data)[ki];
2125
2126      //Quantities at neighbour on nearest face
2127      n = ((long *) neighbours -> data)[ki];
2128      if (n < 0) {
2129        m = -n-1; //Convert negative flag to index
2130        qr[0] = ((double *) stage_boundary_values -> data)[m];
2131        qr[1] = ((double *) xmom_boundary_values -> data)[m];
2132        qr[2] = ((double *) ymom_boundary_values -> data)[m];
2133        zr = zl; //Extend bed elevation to boundary
2134      } else {
2135        m = ((long *) neighbour_edges -> data)[ki];
2136        nm = n*3+m;
2137        qr[0] = ((double *) stage_edge_values -> data)[nm];
2138        qr[1] = ((double *) xmom_edge_values -> data)[nm];
2139        qr[2] = ((double *) ymom_edge_values -> data)[nm];
2140        zr =    ((double *) bed_edge_values -> data)[nm];
2141      }
2142      // Outward pointing normal vector
2143      // normal = domain.normals[k, 2*i:2*i+2]
2144      ki2 = 2*ki; //k*6 + i*2
2145      normal[0] = ((double *) normals -> data)[ki2];
2146      normal[1] = ((double *) normals -> data)[ki2+1];
2147      //Edge flux computation
2148      flux_function_kinetic(ql, qr, zl, zr,
2149                    normal[0], normal[1],
2150                    epsilon, H0, g,
2151                    edgeflux, &max_speed);
2152      //update triangle k
2153      ((long *) already_computed_flux->data)[ki]=call;
2154      ((double *) stage_explicit_update -> data)[k] -= edgeflux[0]*((double *) edgelengths -> data)[ki];
2155      ((double *) xmom_explicit_update -> data)[k] -= edgeflux[1]*((double *) edgelengths -> data)[ki];
2156      ((double *) ymom_explicit_update -> data)[k] -= edgeflux[2]*((double *) edgelengths -> data)[ki];
2157      //update the neighbour n
2158      if (n>=0){
2159        ((long *) already_computed_flux->data)[nm]=call;
2160        ((double *) stage_explicit_update -> data)[n] += edgeflux[0]*((double *) edgelengths -> data)[nm];
2161        ((double *) xmom_explicit_update -> data)[n] += edgeflux[1]*((double *) edgelengths -> data)[nm];
2162        ((double *) ymom_explicit_update -> data)[n] += edgeflux[2]*((double *) edgelengths -> data)[nm];
2163      }
2164      ///for (j=0; j<3; j++) {
2165        ///flux[j] -= edgeflux[j]*((double *) edgelengths -> data)[ki];
2166        ///}
2167        //Update timestep
2168        //timestep = min(timestep, domain.radii[k]/max_speed)
2169        //FIXME: SR Add parameter for CFL condition
2170    if ( ((long *) tri_full_flag -> data)[k] == 1) {
2171            if (max_speed > epsilon) {
2172                timestep = min(timestep, ((double *) radii -> data)[k]/max_speed);
2173                //maxspeed in flux_function is calculated as max(|u+a|,|u-a|)
2174                if (n>=0)
2175                    timestep = min(timestep, ((double *) radii -> data)[n]/max_speed);
2176            }
2177    }
2178    } // end for i
2179    //Normalise by area and store for when all conserved
2180    //quantities get updated
2181    ((double *) stage_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
2182    ((double *) xmom_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
2183    ((double *) ymom_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
2184  } //end for k
2185  return Py_BuildValue("d", timestep);
2186}
2187
2188PyObject *protect(PyObject *self, PyObject *args) {
2189  //
2190  //    protect(minimum_allowed_height, maximum_allowed_speed, wc, zc, xmomc, ymomc)
2191
2192
2193  PyArrayObject
2194  *wc,            //Stage at centroids
2195  *zc,            //Elevation at centroids
2196  *xmomc,         //Momentums at centroids
2197  *ymomc;
2198
2199
2200  int N;
2201  double minimum_allowed_height, maximum_allowed_speed, epsilon;
2202
2203  // Convert Python arguments to C
2204  if (!PyArg_ParseTuple(args, "dddOOOO",
2205                        &minimum_allowed_height,
2206                        &maximum_allowed_speed,
2207                        &epsilon,
2208                        &wc, &zc, &xmomc, &ymomc)) {
2209    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: protect could not parse input arguments");
2210    return NULL;
2211  } 
2212
2213  N = wc -> dimensions[0];
2214
2215  _protect(N,
2216           minimum_allowed_height,
2217           maximum_allowed_speed,
2218           epsilon,
2219           (double*) wc -> data,
2220           (double*) zc -> data,
2221           (double*) xmomc -> data,
2222           (double*) ymomc -> data);
2223
2224  return Py_BuildValue("");
2225}
2226
2227
2228
2229PyObject *balance_deep_and_shallow(PyObject *self, PyObject *args) {
[4769]2230  // Compute linear combination between stage as computed by
2231  // gradient-limiters limiting using w, and stage computed by
2232  // gradient-limiters limiting using h (h-limiter).
2233  // The former takes precedence when heights are large compared to the
2234  // bed slope while the latter takes precedence when heights are
2235  // relatively small.  Anything in between is computed as a balanced
2236  // linear combination in order to avoid numerical disturbances which
2237  // would otherwise appear as a result of hard switching between
2238  // modes.
[4376]2239  //
[4733]2240  //    balance_deep_and_shallow(beta_h, wc, zc, wv, zv,
[4376]2241  //                             xmomc, ymomc, xmomv, ymomv)
2242
2243
2244  PyArrayObject
2245    *wc,            //Stage at centroids
2246    *zc,            //Elevation at centroids
2247    *wv,            //Stage at vertices
2248    *zv,            //Elevation at vertices
2249    *hvbar,         //h-Limited depths at vertices
2250    *xmomc,         //Momentums at centroids and vertices
2251    *ymomc,
2252    *xmomv,
2253    *ymomv;
2254 
2255  PyObject *domain, *Tmp;
2256   
2257  double alpha_balance = 2.0;
[4733]2258  double H0, beta_h;
[4376]2259
[4631]2260  int N, tight_slope_limiters; //, err;
[4376]2261
2262  // Convert Python arguments to C
[4733]2263  if (!PyArg_ParseTuple(args, "OdOOOOOOOOO",
[4376]2264                        &domain,
[4733]2265                        &beta_h,
2266                        &wc, &zc, 
2267                        &wv, &zv, &hvbar,
[4376]2268                        &xmomc, &ymomc, &xmomv, &ymomv)) {
[4733]2269    PyErr_SetString(PyExc_RuntimeError, 
2270                    "shallow_water_ext.c: balance_deep_and_shallow could not parse input arguments");
[4376]2271    return NULL;
2272  } 
2273         
[4730]2274         
2275  // FIXME (Ole): I tested this without GetAttrString and got time down
2276  // marginally from 4.0s to 3.8s. Consider passing everything in
2277  // through ParseTuple and profile.
2278 
[4376]2279  // Pull out parameters
2280  Tmp = PyObject_GetAttrString(domain, "alpha_balance");
2281  if (!Tmp) {
2282    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not obtain object alpha_balance from domain");
2283    return NULL;
2284  } 
2285  alpha_balance = PyFloat_AsDouble(Tmp);
2286  Py_DECREF(Tmp);
2287
2288 
2289  Tmp = PyObject_GetAttrString(domain, "H0");
2290  if (!Tmp) {
2291    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not obtain object H0 from domain");
2292    return NULL;
2293  } 
2294  H0 = PyFloat_AsDouble(Tmp);
2295  Py_DECREF(Tmp);
2296
2297 
[4631]2298  Tmp = PyObject_GetAttrString(domain, "tight_slope_limiters");
[4376]2299  if (!Tmp) {
[4631]2300    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not obtain object tight_slope_limiters from domain");
[4376]2301    return NULL;
2302  } 
[4631]2303  tight_slope_limiters = PyInt_AsLong(Tmp);
[4376]2304  Py_DECREF(Tmp);
2305 
2306 
2307  N = wc -> dimensions[0];
2308  _balance_deep_and_shallow(N,
[4733]2309                            beta_h,
[4376]2310                            (double*) wc -> data,
2311                            (double*) zc -> data,
2312                            (double*) wv -> data,
2313                            (double*) zv -> data,
2314                            (double*) hvbar -> data,
2315                            (double*) xmomc -> data,
2316                            (double*) ymomc -> data,
2317                            (double*) xmomv -> data,
2318                            (double*) ymomv -> data,
2319                            H0,
[4631]2320                            (int) tight_slope_limiters,
[4376]2321                            alpha_balance);
2322
2323
2324  return Py_BuildValue("");
2325}
2326
2327
2328
2329PyObject *h_limiter(PyObject *self, PyObject *args) {
2330
2331  PyObject *domain, *Tmp;
2332  PyArrayObject
2333    *hv, *hc, //Depth at vertices and centroids
2334    *hvbar,   //Limited depth at vertices (return values)
2335    *neighbours;
2336
2337  int k, i, n, N, k3;
2338  int dimensions[2];
[4769]2339  double beta_h; // Safety factor (see config.py)
[4376]2340  double *hmin, *hmax, hn;
2341
2342  // Convert Python arguments to C
2343  if (!PyArg_ParseTuple(args, "OOO", &domain, &hc, &hv)) {
2344    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: h_limiter could not parse input arguments");
2345    return NULL;
2346  } 
2347 
2348  neighbours = get_consecutive_array(domain, "neighbours");
2349
2350  //Get safety factor beta_h
2351  Tmp = PyObject_GetAttrString(domain, "beta_h");
2352  if (!Tmp) {
2353    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: h_limiter could not obtain object beta_h from domain");
2354    return NULL;
2355  } 
2356  beta_h = PyFloat_AsDouble(Tmp);
2357
2358  Py_DECREF(Tmp);
2359
2360  N = hc -> dimensions[0];
2361
2362  //Create hvbar
2363  dimensions[0] = N;
2364  dimensions[1] = 3;
2365  hvbar = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
2366
2367
2368  //Find min and max of this and neighbour's centroid values
2369  hmin = malloc(N * sizeof(double));
2370  hmax = malloc(N * sizeof(double));
2371  for (k=0; k<N; k++) {
2372    k3=k*3;
2373
2374    hmin[k] = ((double*) hc -> data)[k];
2375    hmax[k] = hmin[k];
2376
2377    for (i=0; i<3; i++) {
2378      n = ((long*) neighbours -> data)[k3+i];
2379
[4769]2380      // Initialise hvbar with values from hv
[4376]2381      ((double*) hvbar -> data)[k3+i] = ((double*) hv -> data)[k3+i];
2382
2383      if (n >= 0) {
[4769]2384        hn = ((double*) hc -> data)[n]; // Neighbour's centroid value
[4376]2385
2386        hmin[k] = min(hmin[k], hn);
2387        hmax[k] = max(hmax[k], hn);
2388      }
2389    }
2390  }
2391
2392  // Call underlying standard routine
[4883]2393  _limit_old(N, beta_h, (double*) hc -> data, (double*) hvbar -> data, hmin, hmax);
[4376]2394
[4769]2395  // // //Py_DECREF(domain); //FIXME: Necessary?
[4376]2396  free(hmin);
2397  free(hmax);
2398
[4769]2399  // Return result using PyArray to avoid memory leak
[4376]2400  return PyArray_Return(hvbar);
2401}
2402
2403PyObject *h_limiter_sw(PyObject *self, PyObject *args) {
2404  //a faster version of h_limiter above
2405  PyObject *domain, *Tmp;
2406  PyArrayObject
2407    *hv, *hc, //Depth at vertices and centroids
2408    *hvbar,   //Limited depth at vertices (return values)
2409    *neighbours;
2410
2411  int k, i, N, k3,k0,k1,k2;
2412  int dimensions[2];
2413  double beta_h; //Safety factor (see config.py)
2414  double hmin, hmax, dh[3];
2415// Convert Python arguments to C
2416  if (!PyArg_ParseTuple(args, "OOO", &domain, &hc, &hv)) {
2417    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: h_limiter_sw could not parse input arguments");
2418    return NULL;
2419  } 
2420  neighbours = get_consecutive_array(domain, "neighbours");
2421
2422  //Get safety factor beta_h
2423  Tmp = PyObject_GetAttrString(domain, "beta_h");
2424  if (!Tmp) {
2425    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: h_limiter_sw could not obtain object beta_h from domain");
2426    return NULL;
2427  } 
2428  beta_h = PyFloat_AsDouble(Tmp);
2429
2430  Py_DECREF(Tmp);
2431
2432  N = hc -> dimensions[0];
2433
2434  //Create hvbar
2435  dimensions[0] = N;
2436  dimensions[1] = 3;
2437  hvbar = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
2438  for (k=0;k<N;k++){
2439    k3=k*3;
2440    //get the ids of the neighbours
2441    k0 = ((long*) neighbours -> data)[k3];
2442    k1 = ((long*) neighbours -> data)[k3+1];
2443    k2 = ((long*) neighbours -> data)[k3+2];
2444    //set hvbar provisionally
2445    for (i=0;i<3;i++){
2446      ((double*) hvbar -> data)[k3+i] = ((double*) hv -> data)[k3+i];
2447      dh[i]=((double*) hvbar -> data)[k3+i]-((double*) hc -> data)[k];
2448    }
2449    hmin=((double*) hc -> data)[k];
2450    hmax=hmin;
2451    if (k0>=0){
2452      hmin=min(hmin,((double*) hc -> data)[k0]);
2453      hmax=max(hmax,((double*) hc -> data)[k0]);
2454    }
2455    if (k1>=0){
2456      hmin=min(hmin,((double*) hc -> data)[k1]);
2457      hmax=max(hmax,((double*) hc -> data)[k1]);
2458    }
2459    if (k2>=0){
2460      hmin=min(hmin,((double*) hc -> data)[k2]);
2461      hmax=max(hmax,((double*) hc -> data)[k2]);
2462    }
2463    hmin-=((double*) hc -> data)[k];
2464    hmax-=((double*) hc -> data)[k];
2465    limit_gradient(dh,hmin,hmax,beta_h);
2466    for (i=0;i<3;i++)
2467      ((double*) hvbar -> data)[k3+i] = ((double*) hc -> data)[k]+dh[i];
2468  }
2469  return PyArray_Return(hvbar);
2470}
2471
2472PyObject *assign_windfield_values(PyObject *self, PyObject *args) {
2473  //
2474  //      assign_windfield_values(xmom_update, ymom_update,
2475  //                              s_vec, phi_vec, self.const)
2476
2477
2478
2479  PyArrayObject   //(one element per triangle)
2480  *s_vec,         //Speeds
2481  *phi_vec,       //Bearings
2482  *xmom_update,   //Momentum updates
2483  *ymom_update;
2484
2485
2486  int N;
2487  double cw;
2488
2489  // Convert Python arguments to C
2490  if (!PyArg_ParseTuple(args, "OOOOd",
2491                        &xmom_update,
2492                        &ymom_update,
2493                        &s_vec, &phi_vec,
2494                        &cw)) {
2495    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: assign_windfield_values could not parse input arguments");
2496    return NULL;
2497  } 
2498                       
2499
2500  N = xmom_update -> dimensions[0];
2501
2502  _assign_wind_field_values(N,
2503           (double*) xmom_update -> data,
2504           (double*) ymom_update -> data,
2505           (double*) s_vec -> data,
2506           (double*) phi_vec -> data,
2507           cw);
2508
2509  return Py_BuildValue("");
2510}
2511
2512
2513
2514
[4769]2515//-------------------------------
[4376]2516// Method table for python module
[4769]2517//-------------------------------
[4376]2518static struct PyMethodDef MethodTable[] = {
2519  /* The cast of the function is necessary since PyCFunction values
2520   * only take two PyObject* parameters, and rotate() takes
2521   * three.
2522   */
2523
2524  {"rotate", (PyCFunction)rotate, METH_VARARGS | METH_KEYWORDS, "Print out"},
2525  {"extrapolate_second_order_sw", extrapolate_second_order_sw, METH_VARARGS, "Print out"},
2526  {"compute_fluxes_ext_central", compute_fluxes_ext_central, METH_VARARGS, "Print out"},
2527  {"compute_fluxes_ext_kinetic", compute_fluxes_ext_kinetic, METH_VARARGS, "Print out"},
2528  {"gravity", gravity, METH_VARARGS, "Print out"},
2529  {"manning_friction", manning_friction, METH_VARARGS, "Print out"},
[4769]2530  {"flux_function_central", flux_function_central, METH_VARARGS, "Print out"}, 
[4376]2531  {"balance_deep_and_shallow", balance_deep_and_shallow,
2532   METH_VARARGS, "Print out"},
2533  {"h_limiter", h_limiter,
2534   METH_VARARGS, "Print out"},
2535  {"h_limiter_sw", h_limiter_sw,
2536   METH_VARARGS, "Print out"},
2537  {"protect", protect, METH_VARARGS | METH_KEYWORDS, "Print out"},
2538  {"assign_windfield_values", assign_windfield_values,
2539   METH_VARARGS | METH_KEYWORDS, "Print out"},
2540  {NULL, NULL}
2541};
2542
2543// Module initialisation
2544void initshallow_water_ext(void){
2545  Py_InitModule("shallow_water_ext", MethodTable);
2546
[4733]2547  import_array(); // Necessary for handling of NumPY structures
[4376]2548}
Note: See TracBrowser for help on using the repository browser.