# source:anuga_core/source/anuga/shallow_water/shallow_water_ext.c@7034

Last change on this file since 7034 was 7034, checked in by ole, 14 years ago

Comments on the fast sqrt approximations. They don't give speedup on AMD 64, so have been left alone until further notice.
Timings included in the code.

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