# source:trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain_ext.c@8446

Last change on this file since 8446 was 8446, checked in by davies, 11 years ago

File size: 66.7 KB
Line
1// Python - C extension module for shallow_water.py
2//
3// To compile (Python2.6):
4//  gcc -c swb2_domain_ext.c -I/usr/include/python2.6 -o domain_ext.o -Wall -O
5//  gcc -shared swb2_domain_ext.o  -o swb2_domain_ext.so
6//
7// or use python compile.py
8//
9// See the module swb_domain.py for more documentation on
10// how to use this module
11//
12//
13// Stephen Roberts, ANU 2009
14// Ole Nielsen, GA 2004
15// Gareth Davies, GA 2011
16
17#include "Python.h"
18#include "numpy/arrayobject.h"
19#include "math.h"
20#include <stdio.h>
21//#include "numpy_shim.h"
22
23// Shared code snippets
24#include "util_ext.h"
25
26
27const double pi = 3.14159265358979;
28
29// Computational function for rotation
30int _rotate(double *q, double n1, double n2) {
31  /*Rotate the last  2 coordinates of 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 2x1 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
42  // Shorthands
43  q1 = q[1];  // x coordinate
44  q2 = q[2];  // y coordinate
45
46  // Rotate
47  q[1] =  n1*q1 + n2*q2;
48  q[2] = -n2*q1 + n1*q2;
49
50  return 0;
51}
52
53// Function to obtain speed from momentum and depth.
54// This is used by flux functions
55// Input parameters uh and h may be modified by this function.
56// Tried to inline, but no speedup was achieved 27th May 2009 (Ole)
57//static inline double _compute_speed(double *uh,
58double _compute_speed(double *uh,
59                      double *h,
60                      double epsilon,
61                      double h0,
62                      double limiting_threshold) {
63
64  double u;
65
66  if (*h < limiting_threshold) {
67    // Apply limiting of speeds according to the ANUGA manual
68    if (*h < epsilon) {
69      //*h = max(0.0,*h);  // Could have been negative
70      *h = 0.0;  // Could have been negative
71      u = 0.0;
72    } else {
73      u = *uh/(*h + h0/ *h);
74    }
75
76
77    // Adjust momentum to be consistent with speed
78    *uh = u * *h;
79  } else {
80    // We are in deep water - no need for limiting
81    u = *uh/ *h;
82  }
83
84  return u;
85}
86
87// Innermost flux function (using stage w=z+h)
88int _flux_function_central(double *q_left, double *q_right,
89                           double z_left, double z_right,
90                           double n1, double n2,
91                           double epsilon,
92                           double h0,
93                           double limiting_threshold,
94                           double g,
95                           double *edgeflux, double *max_speed,
96                           double *pressure_flux)
97{
98
99  /*Compute fluxes between volumes for the shallow water wave equation
100    cast in terms of the 'stage', w = h+z using
101    the 'central scheme' as described in
102
103    Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
104    Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
105    Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.
106
107    The implemented formula is given in equation (3.15) on page 714
108  */
109
110  int i;
111
112  double w_left, h_left, uh_left, vh_left, u_left;
113  double w_right, h_right, uh_right, vh_right, u_right;
114  double s_min, s_max, soundspeed_left, soundspeed_right;
115  double denom, inverse_denominator, z;
116  double uint, t1, t2, t3;
117  // Workspace (allocate once, use many)
118  static double q_left_rotated[3], q_right_rotated[3], flux_right[3], flux_left[3];
119
120  // Copy conserved quantities to protect from modification
121  q_left_rotated[0] = q_left[0];
122  q_right_rotated[0] = q_right[0];
123  q_left_rotated[1] = q_left[1];
124  q_right_rotated[1] = q_right[1];
125  q_left_rotated[2] = q_left[2];
126  q_right_rotated[2] = q_right[2];
127
128  // Align x- and y-momentum with x-axis
129  _rotate(q_left_rotated, n1, n2);
130  _rotate(q_right_rotated, n1, n2);
131
132  z = 0.5*(z_left + z_right); // Average elevation values.
133                              // Even though this will nominally allow
134                              // for discontinuities in the elevation data,
135                              // there is currently no numerical support for
136                              // this so results may be strange near
137                              // jumps in the bed.
138
139  // Compute speeds in x-direction
140  w_left = q_left_rotated[0];
141  h_left = w_left - z;
142  uh_left = q_left_rotated[1];
143  u_left = _compute_speed(&uh_left, &h_left,
144                          epsilon, h0, limiting_threshold);
145
146  w_right = q_right_rotated[0];
147  h_right = w_right - z;
148  uh_right = q_right_rotated[1];
149  u_right = _compute_speed(&uh_right, &h_right,
150                           epsilon, h0, limiting_threshold);
151
152  // Momentum in y-direction
153  vh_left  = q_left_rotated[2];
154  vh_right = q_right_rotated[2];
155
156  // Limit y-momentum if necessary
157  // Leaving this out, improves speed significantly (Ole 27/5/2009)
158  // All validation tests pass, so do we really need it anymore?
159  _compute_speed(&vh_left, &h_left,
160                 epsilon, h0, limiting_threshold);
161  _compute_speed(&vh_right, &h_right,
162                 epsilon, h0, limiting_threshold);
163
164  // Maximal and minimal wave speeds
165  soundspeed_left  = sqrt(g*h_left);
166  soundspeed_right = sqrt(g*h_right);
167
168  // Code to use fast square root optimisation if desired.
169  // Timings on AMD 64 for the Okushiri profile gave the following timings
170  //
171  // SQRT           Total    Flux
172  //=============================
173  //
174  // Ref            405s     152s
175  // Fast (dbl)     453s     173s
176  // Fast (sng)     437s     171s
177  //
178  // Consequently, there is currently (14/5/2009) no reason to use this
179  // approximation.
180
181  //soundspeed_left  = fast_squareroot_approximation(g*h_left);
182  //soundspeed_right = fast_squareroot_approximation(g*h_right);
183
184  s_max = max(u_left + soundspeed_left, u_right + soundspeed_right);
185  if (s_max < 0.0)
186  {
187    s_max = 0.0;
188  }
189
190  s_min = min(u_left - soundspeed_left, u_right - soundspeed_right);
191  if (s_min > 0.0)
192  {
193    s_min = 0.0;
194  }
195
196  // Flux formulas
197  flux_left[0] = u_left*h_left;
198  flux_left[1] = u_left*uh_left + 0.5*g*h_left*h_left;
199  flux_left[2] = u_left*vh_left;
200
201  flux_right[0] = u_right*h_right;
202  flux_right[1] = u_right*uh_right + 0.5*g*h_right*h_right;
203  flux_right[2] = u_right*vh_right;
204
205  // Flux computation
206  denom = s_max - s_min;
207  if (denom < epsilon)
208  { // FIXME (Ole): Try using h0 here
209    memset(edgeflux, 0, 3*sizeof(double));
210    //for(i=0;i<3;i++){
211    //    edgeflux[i] = _minmod(flux_left[i], flux_right[i]);
212    //}
213    *max_speed = 0.0;
214  }
215  else
216  {
217    inverse_denominator = 1.0/denom;
218    for (i = 0; i < 3; i++)
219    {
220      // Adjustment to the scheme by Kurganov and Lin 2007 Communications in Computational
221      // Physics 2:141-163
222      //uint = (s_max*q_right_rotated[i] - s_min*q_left_rotated[i] - (flux_right[i] - flux_left[i]))*inverse_denominator;
223      //t1 = (q_right_rotated[i] - uint);
224      //t2 = (-q_left_rotated[i] + uint);
225      //t3 = _minmod(t1,t2);
226      t3 = 0.0;
227      edgeflux[i] = s_max*flux_left[i] - s_min*flux_right[i];
228      //if(i<2) {
229      edgeflux[i] += s_max*s_min*(q_right_rotated[i] - q_left_rotated[i] - t3);
230      //}else{
231      //  edgeflux[i] += 0.5*s_max*s_min*(q_right_rotated[i] - q_left_rotated[i] - t3);
232      //}
233      //t1 = s_max*s_min*(q_right_rotated[i] - q_left_rotated[i] - t3);
234      //if(fabs(edgeflux[i])>fabs(t1)){
235      //      edgeflux[i]+=t1;
236      //}else{
237      //      edgeflux[i]+=fabs(edgeflux[i])*copysign(1.0,t1) ;
238      //}
239
240      edgeflux[i] *= inverse_denominator;
241    }
242
243    *pressure_flux = 0.5*g*( s_max*h_left*h_left -s_min*h_right*h_right)*inverse_denominator;
244    //if(edgeflux[0]<0.0){
245    //    //edgeflux[2] = edgeflux[0]*0.5*(vh_right/h_right + vh_left/h_left);
246    //    edgeflux[2] = edgeflux[0]*vh_right/h_right;
247    //}else{
248    //    edgeflux[2] = edgeflux[0]*vh_left/h_left;
249    //}
250
251    // Maximal wavespeed
252    *max_speed = max(fabs(s_max), fabs(s_min));
253
254    // Rotate back
255    _rotate(edgeflux, n1, -n2);
256  }
257
258  return 0;
259}
260
261
262// Computational function for flux computation
263double _compute_fluxes_central(int number_of_elements,
264        double timestep,
265        double epsilon,
266        double H0,
267        double g,
268        long* neighbours,
269        long* neighbour_edges,
270        double* normals,
271        double* edgelengths,
273        double* areas,
274        long* tri_full_flag,
275        double* stage_edge_values,
276        double* xmom_edge_values,
277        double* ymom_edge_values,
278        double* bed_edge_values,
279        double* stage_boundary_values,
280        double* xmom_boundary_values,
281        double* ymom_boundary_values,
282        long*   boundary_flux_type,
283        double* stage_explicit_update,
284        double* xmom_explicit_update,
285        double* ymom_explicit_update,
287        double* max_speed_array,
288        int optimise_dry_cells,
289        double* stage_centroid_values,
290        double* bed_centroid_values,
291        double* bed_vertex_values) {
292    // Local variables
293    double max_speed, length, inv_area, zl, zr;
294    double h0 = H0*H0; // This ensures a good balance when h approaches H0.
295
296    double limiting_threshold = 10 * H0; // Avoid applying limiter below this
297    // threshold for performance reasons.
298    // See ANUGA manual under flux limiting
299    int k, i, m, n,j;
300    int ki, nm = 0, ki2; // Index shorthands
301
302    // Workspace (making them static actually made function slightly slower (Ole))
303    double ql[3], qr[3], edgeflux[3]; // Work array for summing up fluxes
304    double stage_edges[3];//Work array
305    double bedslope_work;
306    int neighbours_wet[3];//Work array
307    int useint;
308    double stage_edge_lim, scale_factor_shallow, bedtop, bedbot, pressure_flux, hc, hc_n;
309    //double min_bed_edgevalue[number_of_elements]; // FIXME: Set memory explicitly
310    static long call = 1; // Static local variable flagging already computed flux
311
312    double *min_bed_edgevalue;
313
314    min_bed_edgevalue = malloc(number_of_elements*sizeof(double));
315    // Start computation
316    call++; // Flag 'id' of flux calculation for this timestep
317
318    // Set explicit_update to zero for all conserved_quantities.
319    // This assumes compute_fluxes called before forcing terms
320    memset((char*) stage_explicit_update, 0, number_of_elements * sizeof (double));
321    memset((char*) xmom_explicit_update, 0, number_of_elements * sizeof (double));
322    memset((char*) ymom_explicit_update, 0, number_of_elements * sizeof (double));
323
324
325    // Compute minimum bed edge value on each triangle
326    for (k = 0; k < number_of_elements; k++){
327        min_bed_edgevalue[k] = min(bed_edge_values[3*k],
328                                   min(bed_edge_values[3*k+1], bed_edge_values[3*k+2]));
329        //min_bed_edgevalue[k] = bed_centroid_values[k];
330
331    }
332
333
334    // For all triangles
335    for (k = 0; k < number_of_elements; k++) {
336        // Loop through neighbours and compute edge flux for each
337        for (i = 0; i < 3; i++) {
338            ki = k * 3 + i; // Linear index to edge i of triangle k
339
340            if (already_computed_flux[ki] == call) {
341                // We've already computed the flux across this edge
342                continue;
343            }
344
345            // Get left hand side values from triangle k, edge i
346            ql[0] = stage_edge_values[ki];
347            ql[1] = xmom_edge_values[ki];
348            ql[2] = ymom_edge_values[ki];
349            zl = bed_edge_values[ki];
350            hc = max(stage_centroid_values[k] - bed_centroid_values[k],0.0);
351            // Get right hand side values either from neighbouring triangle
352            // or from boundary array (Quantities at neighbour on nearest face).
353            n = neighbours[ki];
354            hc_n = hc;
355            if (n < 0) {
356                // Neighbour is a boundary condition
357                m = -n - 1; // Convert negative flag to boundary index
358
359                qr[0] = stage_boundary_values[m];
360                qr[1] = xmom_boundary_values[m];
361                qr[2] = ymom_boundary_values[m];
362                zr = zl; // Extend bed elevation to boundary
363            }
364            else {
365                // Neighbour is a real triangle
366                hc_n = max(stage_centroid_values[n] - bed_centroid_values[n],0.0);
367                m = neighbour_edges[ki];
368                nm = n * 3 + m; // Linear index (triangle n, edge m)
369
370                qr[0] = stage_edge_values[nm];
371                qr[1] = xmom_edge_values[nm];
372                qr[2] = ymom_edge_values[nm];
373                zr = bed_edge_values[nm];
374            }
375
376            // Now we have values for this edge - both from left and right side.
377
378            if (optimise_dry_cells) {
379                // Check if flux calculation is necessary across this edge
380                // This check will exclude dry cells.
381                // This will also optimise cases where zl != zr as
382                // long as both are dry
383
384                if ((ql[0] - zl) < epsilon &&
385                        (qr[0] - zr) < epsilon) {
386                    // Cell boundary is dry
387
388                    already_computed_flux[ki] = call; // #k Done
389                    if (n >= 0) {
390                        already_computed_flux[nm] = call; // #n Done
391                    }
392
393                    max_speed = 0.0;
394                    continue;
395                }
396            }
397
398
399            if (fabs(zl-zr)>1.0e-10) {
400                report_python_error(AT,"Discontinuous Elevation");
401                return 0.0;
402            }
403
404            // If both centroids are dry, then the cells should not exchange flux
405            // NOTE: This also means no bedslope term -- which is arguably
406            // appropriate, since the depth is zero at the cell centroid
407            if(( hc == 0.0)&(hc_n == 0.0)){
408                already_computed_flux[ki] = call; // #k Done
409                already_computed_flux[nm] = call; // #n Done
410                max_speed = 0.0;
411                continue;
412            }
413
414            //If one centroid is dry, then extrapolate its edge values from the neighbouring centroid,
415            // unless the local centroid value is smaller
416            if(n>=0){
417                if(hc==0.0){
418                    //ql[0]=max(min(qr[0],stage_centroid_values[k]),zl);
419                    //ql[0]=max(min(qr[0],0.5*(stage_centroid_values[k]+stage_centroid_values[n])),zl);
420                    ql[0]=qr[0]; // max(min(qr[0],stage_centroid_values[k]),zl);
421                }
422                if(hc_n==0.0){
423                    //qr[0]=max(min(ql[0],stage_centroid_values[n]),zr);
424                    //qr[0]=max(min(ql[0],0.5*(stage_centroid_values[n]+stage_centroid_values[k])),zr);
425                    qr[0]=ql[0];
426                }
427            }else{
428                // Treat the boundary case
429                if((hc==0.0)){
430                    ql[0]=max(min(qr[0],stage_centroid_values[k]),zl);
431                    //ql[0]=max(min(qr[0],ql[0]),zl);
432                }
433            }
434
435            // Outward pointing normal vector (domain.normals[k, 2*i:2*i+2])
436            ki2 = 2 * ki; //k*6 + i*2
437
438            // Edge flux computation (triangle k, edge i)
439            _flux_function_central(ql, qr, zl, zr,
440                    normals[ki2], normals[ki2 + 1],
441                    epsilon, h0, limiting_threshold, g,
442                    edgeflux, &max_speed, &pressure_flux);
443
444            // Prevent outflow from 'seriously' dry cells
445            if(stage_centroid_values[k]<min_bed_edgevalue[k]){
446                if(edgeflux[0]>0.0){
447                    edgeflux[0]=0.0;
448                    edgeflux[1]=0.0;
449                    edgeflux[2]=0.0;
450                }
451            }
452            if(n>=0){
453                if(stage_centroid_values[n]<min_bed_edgevalue[n]){
454                    if(edgeflux[0]<0.0){
455                        edgeflux[0]=0.0;
456                        edgeflux[1]=0.0;
457                        edgeflux[2]=0.0;
458                    }
459                }
460            }
461
462            // Multiply edgeflux by edgelength
463            length = edgelengths[ki];
464            edgeflux[0] *= length;
465            edgeflux[1] *= length;
466            edgeflux[2] *= length;
467
468            // GET RID OF THIS: EXPERIMENTAL
469            //if(n<0){
470            //    if(boundary_flux_type[m]>0){
471            //        // IMPOSE BOUNDARY CONDITIONS DIRECTLY ON THE FLUX
472            //        if(boundary_flux_type[m]==1){
473            //            // Zero_mass_flux_transmissive_momentum_flux boundary, or
474            //            // zero_mass_flux_zero_momentum_boundary
475            //            // Just need to set the mass flux to zero, as the
476            //            // transmissive momentum is ensured by the values of
477            //            // qr[0], qr[1], qr[2]
478            //            printf("Warning: WEIRD EDGEFLUX THING IS ACTING");
479            //            edgeflux[0] = 0.0;
480            //        }
481            //    }
482            //}
483
484            // Update triangle k with flux from edge i
485            stage_explicit_update[k] -= edgeflux[0];
486            xmom_explicit_update[k] -= edgeflux[1];
487            ymom_explicit_update[k] -= edgeflux[2];
488
489            // Compute bed slope term
490            if(hc>-9999.0){
491                //Bedslope approx 1:
492                //bedslope_work = g*hc*(ql[0])*length-0.5*g*max(ql[0]-zl,0.)*(ql[0]-zl)*length;
493                //
494                // Bedslope approx 2
495                //stage_edge_lim = ql[0]; // Limit this to be between a constant stage and constant depth extrapolation
496                //if(stage_edge_lim > max(stage_centroid_values[k], zl +hc)){
497                //    stage_edge_lim = max(stage_centroid_values[k], zl+hc);
498                //}
499                //if(stage_edge_lim < min(stage_centroid_values[k], zl +hc)){
500                //    stage_edge_lim = min(stage_centroid_values[k], zl+hc);
501                //}
502                //bedslope_work = g*hc*(stage_edge_lim)*length-0.5*g*max(stage_edge_lim-zl,0.)*(stage_edge_lim-zl)*length;
503
504                // Bedslope approx 3
505                bedslope_work = -0.5*g*max(stage_centroid_values[k]-zl,0.)*(stage_centroid_values[k]-zl)*length;
506                //
507                xmom_explicit_update[k] -= normals[ki2]*bedslope_work;
508                ymom_explicit_update[k] -= normals[ki2+1]*bedslope_work;
509            }else{
510                // Treat nearly dry cells
511                bedslope_work =-0.5*g*max(ql[0]-zl,0.)*(ql[0]-zl)*length; //
512                //
513                //bedslope_work = -pressure_flux*length;
514                xmom_explicit_update[k] -= normals[ki2]*bedslope_work;
515                ymom_explicit_update[k] -= normals[ki2+1]*bedslope_work;
516
517                //xmom_explicit_update[k] = 0.0;
518                //ymom_explicit_update[k] = 0.0;
519
520            }
521
522            already_computed_flux[ki] = call; // #k Done
523
524
525            // Update neighbour n with same flux but reversed sign
526            if (n >= 0) {
527                stage_explicit_update[n] += edgeflux[0];
528                xmom_explicit_update[n] += edgeflux[1];
529                ymom_explicit_update[n] += edgeflux[2];
530                //Add bed slope term here
531                if(hc_n>-9999.0){
532                //if(stage_centroid_values[n] > max_bed_edgevalue[n]){
533                    // Bedslope approx 1:
534                    //bedslope_work = g*hc_n*(qr[0])*length-0.5*g*max(qr[0]-zr,0.)*(qr[0]-zr)*length;
535                    //
536                    // Bedslope approx 2:
537                    //stage_edge_lim = qr[0];
538                    //if(stage_edge_lim > max(stage_centroid_values[n], zr +hc_n)){
539                    //    stage_edge_lim = max(stage_centroid_values[n], zr+hc_n);
540                    //}
541                    //if(stage_edge_lim < min(stage_centroid_values[n], zr +hc_n)){
542                    //    stage_edge_lim = min(stage_centroid_values[n], zr+hc_n);
543                    //}
544                    //bedslope_work = g*hc_n*(stage_edge_lim)*length-0.5*g*max(stage_edge_lim-zr,0.)*(stage_edge_lim-zr)*length;
545                    //
546                    // Bedslope approx 3
547                    bedslope_work = -0.5*g*max(stage_centroid_values[n]-zr,0.)*(stage_centroid_values[n]-zr)*length;
548                    //
549                    xmom_explicit_update[n] += normals[ki2]*bedslope_work;
550                    ymom_explicit_update[n] += normals[ki2+1]*bedslope_work;
551                }else{
552                    // Treat nearly dry cells
553                    bedslope_work = -0.5*g*max(qr[0]-zr,0.)*(qr[0]-zr)*length; //-pressure_flux*length; //-0.5*g*max(qr[0]-zr,0.)*(qr[0]-zr)*length;
554                    //bedslope_work = -pressure_flux*length;
555                    xmom_explicit_update[n] += normals[ki2]*bedslope_work;
556                    ymom_explicit_update[n] += normals[ki2+1]*bedslope_work;
557
558                    //xmom_explicit_update[n] = 0.0;
559                    //ymom_explicit_update[n] = 0.0;
560                }
561
562                already_computed_flux[nm] = call; // #n Done
563            }
564
565            // Update timestep based on edge i and possibly neighbour n
566            if (tri_full_flag[k] == 1) {
567                if (max_speed > epsilon) {
568                    // Apply CFL condition for triangles joining this edge (triangle k and triangle n)
569
570                    // CFL for triangle k
571                    timestep = min(timestep, radii[k] / max_speed);
572
573                    if (n >= 0) {
574                        // Apply CFL condition for neigbour n (which is on the ith edge of triangle k)
575                        timestep = min(timestep, radii[n] / max_speed);
576                    }
577
578                    // Ted Rigby's suggested less conservative version
579                    //if(n>=0){
581                    //}else{
582                    //  timestep = min(timestep, radii[k]/max_speed);
583                    //}
584                }
585            }
586
587        } // End edge i (and neighbour n)
588
589
590        // Normalise triangle k by area and store for when all conserved
591        // quantities get updated
592        inv_area = 1.0 / areas[k];
593        stage_explicit_update[k] *= inv_area;
594        xmom_explicit_update[k] *= inv_area;
595        ymom_explicit_update[k] *= inv_area;
596
597
598        // Keep track of maximal speeds
599        max_speed_array[k] = max_speed;
600
601    } // End triangle k
602
603    free(min_bed_edgevalue);
604
605    return timestep;
606}
607
608// Protect against the water elevation falling below the triangle bed
609int _protect(int N,
610         double minimum_allowed_height,
611         double maximum_allowed_speed,
612         double epsilon,
613         double* wc,
614         double* wv,
615         double* zc,
616         double* zv,
617         double* xmomc,
618         double* ymomc) {
619
620  int k;
621  double hc, bmin, bmax;
622  double u, v, reduced_speed;
623
624  // This acts like minimum_allowed height, but scales with the vertical
625  // distance between the bed_centroid_value and the max bed_edge_value of
626  // every triangle.
627  double minimum_relative_height=0.1;
628
629
630  // Protect against inifintesimal and negative heights
631  if (maximum_allowed_speed < epsilon) {
632    for (k=0; k<N; k++) {
633      hc = wc[k] - zc[k];
634      // Definine the maximum bed edge value on triangle k.
635      bmax = 0.5*max((zv[3*k]+zv[3*k+1]),max((zv[3*k+1]+zv[3*k+2]),(zv[3*k+2]+zv[3*k])));
636
637      if (hc < max(minimum_relative_height*(bmax-zc[k]), minimum_allowed_height)) {
638
639        // Set momentum to zero and ensure h is non negative
640        // NOTE: THIS IS IMPORTANT -- WE ARE SETTING MOMENTUM TO ZERO
641        xmomc[k] = 0.0;
642        ymomc[k] = 0.0;
643
644        if (hc <= 0.0){
645             // Definine the minimum bed edge value on triangle k.
646             // Setting = minimum edge value can lead to mass conservation problems
647             //bmin = 0.5*min((zv[3*k]+zv[3*k+1]),min((zv[3*k+1]+zv[3*k+2]),(zv[3*k+2]+zv[3*k])));
648             //bmin =0.5*bmin + 0.5*min(zv[3*k],min(zv[3*k+1],zv[3*k+2]));
649             // Setting = minimum vertex value seems better, but might tend to be less smooth
650             bmin =min(zv[3*k],min(zv[3*k+1],zv[3*k+2]));
651             // Minimum allowed stage = bmin
652             // WARNING: ADDING MASS if wc[k]<bmin
653             if(wc[k]<bmin){
654                 printf("Adding mass to dry cell\n");
655             }
656             wc[k] = max(wc[k], bmin) ;
657
658             // Set vertex values as well. Seems that this shouldn't be
659             // needed. However from memory this is important at the first
660             // time step, for 'dry' areas where the designated stage is
661             // less than the bed centroid value
662             wv[3*k] = max(wv[3*k], bmin);
663             wv[3*k+1] = max(wv[3*k+1], bmin);
664             wv[3*k+2] = max(wv[3*k+2], bmin);
665        }
666      }
667    }
668  } else {
669
670    // Protect against infinitesimal and negative heights
671    for (k=0; k<N; k++) {
672      hc = wc[k] - zc[k];
673
674      if (hc < minimum_allowed_height) {
675
676        //New code: Adjust momentum to guarantee speeds are physical
677        //          ensure h is non negative
678
679        if (hc <= 0.0) {
680            wc[k] = zc[k];
681        xmomc[k] = 0.0;
682        ymomc[k] = 0.0;
683        } else {
684          //Reduce excessive speeds derived from division by small hc
685          //FIXME (Ole): This may be unnecessary with new slope limiters
686          //in effect.
687
688          u = xmomc[k]/hc;
689          if (fabs(u) > maximum_allowed_speed) {
690            reduced_speed = maximum_allowed_speed * u/fabs(u);
691            //printf("Speed (u) has been reduced from %.3f to %.3f\n",
692            //   u, reduced_speed);
693            xmomc[k] = reduced_speed * hc;
694          }
695
696          v = ymomc[k]/hc;
697          if (fabs(v) > maximum_allowed_speed) {
698            reduced_speed = maximum_allowed_speed * v/fabs(v);
699            //printf("Speed (v) has been reduced from %.3f to %.3f\n",
700            //   v, reduced_speed);
701            ymomc[k] = reduced_speed * hc;
702          }
703        }
704      }
705    }
706  }
707  return 0;
708}
709
710int find_qmin_and_qmax(double dq0, double dq1, double dq2,
711               double *qmin, double *qmax){
712  // Considering the centroid of an FV triangle and the vertices of its
713  // auxiliary triangle, find
714  // qmin=min(q)-qc and qmax=max(q)-qc,
715  // where min(q) and max(q) are respectively min and max over the
716  // four values (at the centroid of the FV triangle and the auxiliary
717  // triangle vertices),
718  // and qc is the centroid
719  // dq0=q(vertex0)-q(centroid of FV triangle)
720  // dq1=q(vertex1)-q(vertex0)
721  // dq2=q(vertex2)-q(vertex0)
722
723  // This is a simple implementation
724  *qmax = max(max(dq0, max(dq0+dq1, dq0+dq2)), 0.0) ;
725  *qmin = min(min(dq0, min(dq0+dq1, dq0+dq2)), 0.0) ;
726
727  return 0;
728}
729
730int limit_gradient(double *dqv, double qmin, double qmax, double beta_w){
731  // Given provisional jumps dqv from the FV triangle centroid to its
732  // vertices and jumps qmin (qmax) between the centroid of the FV
733  // triangle and the minimum (maximum) of the values at the centroid of
734  // the FV triangle and the auxiliary triangle vertices,
735  // calculate a multiplicative factor phi by which the provisional
736  // vertex jumps are to be limited
737
738  int i;
739  double r=1000.0, r0=1.0, phi=1.0;
740  static double TINY = 1.0e-100; // to avoid machine accuracy problems.
741  // FIXME: Perhaps use the epsilon used elsewhere.
742
743  // Any provisional jump with magnitude < TINY does not contribute to
744  // the limiting process.
745
746  for (i=0;i<3;i++){
747    if (dqv[i]<-TINY)
748      r0=qmin/dqv[i];
749
750    if (dqv[i]>TINY)
751      r0=qmax/dqv[i];
752
753    r=min(r0,r);
754  }
755
756  phi=min(r*beta_w,1.0);
757  //phi=1.;
758  dqv[0]=dqv[0]*phi;
759  dqv[1]=dqv[1]*phi;
760  dqv[2]=dqv[2]*phi;
761
762  return 0;
763}
764
765//int limit_gradient2(double *dqv, double qmin, double qmax, double beta_w, double r0scale){
766//  // EXPERIMENTAL LIMITER, DOESN'T WORK
767//  // Given provisional jumps dqv from the FV triangle centroid to its
768//  // vertices and jumps qmin (qmax) between the centroid of the FV
769//  // triangle and the minimum (maximum) of the values at the centroid of
770//  // the FV triangle and the auxiliary triangle vertices,
771//  // calculate a multiplicative factor phi by which the provisional
772//  // vertex jumps are to be limited
773//
774//  int i;
775//  double r=1000.0, r0=1.0, phi=1.0;
776//  static double TINY = 1.0e-100; // to avoid machine accuracy problems.
777//  // FIXME: Perhaps use the epsilon used elsewhere.
778//
779//  for (i=0;i<3;i++){
780//    if (dqv[i]<-TINY)
781//      r0=qmin/dqv[i]*r0scale*2.0;
782//
783//    if (dqv[i]>TINY)
784//      r0=qmax/dqv[i]*r0scale*2.0;
785//
786//    r=min(r0,r);
787//  }
788//
789//  phi=min(r*beta_w,1.0);
790//  //phi=1.;
791//  //for (i=0;i<3;i++)
792//  dqv[0]=dqv[0]*phi;
793//  dqv[1]=dqv[1]*phi;
794//  dqv[2]=dqv[2]*phi;
795//
796//  return 0;
797//}
798
799// Computational routine
800int _extrapolate_second_order_edge_sw(int number_of_elements,
801                                 double epsilon,
802                                 double minimum_allowed_height,
803                                 double beta_w,
804                                 double beta_w_dry,
805                                 double beta_uh,
806                                 double beta_uh_dry,
807                                 double beta_vh,
808                                 double beta_vh_dry,
809                                 long* surrogate_neighbours,
810                                 long* number_of_boundaries,
811                                 double* centroid_coordinates,
812                                 double* stage_centroid_values,
813                                 double* xmom_centroid_values,
814                                 double* ymom_centroid_values,
815                                 double* elevation_centroid_values,
816                                 double* edge_coordinates,
817                                 double* stage_edge_values,
818                                 double* xmom_edge_values,
819                                 double* ymom_edge_values,
820                                 double* elevation_edge_values,
821                                 double* stage_vertex_values,
822                                 double* xmom_vertex_values,
823                                 double* ymom_vertex_values,
824                                 double* elevation_vertex_values,
825                                 int optimise_dry_cells,
826                                 int extrapolate_velocity_second_order) {
827
828  // Local variables
829  double a, b; // Gradient vector used to calculate edge values from centroids
830  int k, k0, k1, k2, k3, k6, coord_index, i, ii;
831  double x, y, x0, y0, x1, y1, x2, y2, xv0, yv0, xv1, yv1, xv2, yv2; // Vertices of the auxiliary triangle
832  double dx1, dx2, dy1, dy2, dxv0, dxv1, dxv2, dyv0, dyv1, dyv2, dq0, dq1, dq2, area2, inv_area2;
833  double dqv[3], qmin, qmax, hmin, hmax, bedmax, stagemin;
834  double hc, h0, h1, h2, beta_tmp, hfactor;
835  double dk, dv0, dv1, dv2, de[3], demin, dcmax, r0scale;
836
837  double *xmom_centroid_store, *ymom_centroid_store, *stage_centroid_store;
838
839  // Use malloc to avoid putting these variables on the stack, which can cause
840  // segfaults in large model runs
841  xmom_centroid_store = malloc(number_of_elements*sizeof(double));
842  ymom_centroid_store = malloc(number_of_elements*sizeof(double));
843  stage_centroid_store = malloc(number_of_elements*sizeof(double));
844
845  if(extrapolate_velocity_second_order==1){
846      // Replace momentum centroid with velocity centroid to allow velocity
847      // extrapolation This will be changed back at the end of the routine
848      for (k=0; k<number_of_elements; k++){
849
850          dk = max(stage_centroid_values[k]-elevation_centroid_values[k],minimum_allowed_height);
851          xmom_centroid_store[k] = xmom_centroid_values[k];
852          xmom_centroid_values[k] = xmom_centroid_values[k]/dk;
853
854          ymom_centroid_store[k] = ymom_centroid_values[k];
855          ymom_centroid_values[k] = ymom_centroid_values[k]/dk;
856
857      }
858  }
859
860  // Begin extrapolation routine
861  for (k = 0; k < number_of_elements; k++)
862  {
863    k3=k*3;
864    k6=k*6;
865
866    if (number_of_boundaries[k]==3)
867    //if (0==0)
868    {
869      // No neighbours, set gradient on the triangle to zero
870
871      stage_edge_values[k3]   = stage_centroid_values[k];
872      stage_edge_values[k3+1] = stage_centroid_values[k];
873      stage_edge_values[k3+2] = stage_centroid_values[k];
874      xmom_edge_values[k3]    = xmom_centroid_values[k];
875      xmom_edge_values[k3+1]  = xmom_centroid_values[k];
876      xmom_edge_values[k3+2]  = xmom_centroid_values[k];
877      ymom_edge_values[k3]    = ymom_centroid_values[k];
878      ymom_edge_values[k3+1]  = ymom_centroid_values[k];
879      ymom_edge_values[k3+2]  = ymom_centroid_values[k];
880
881      continue;
882    }
883    else
884    {
885      // Triangle k has one or more neighbours.
886      // Get centroid and edge coordinates of the triangle
887
888      // Get the edge coordinates
889      xv0 = edge_coordinates[k6];
890      yv0 = edge_coordinates[k6+1];
891      xv1 = edge_coordinates[k6+2];
892      yv1 = edge_coordinates[k6+3];
893      xv2 = edge_coordinates[k6+4];
894      yv2 = edge_coordinates[k6+5];
895
896      // Get the centroid coordinates
897      coord_index = 2*k;
898      x = centroid_coordinates[coord_index];
899      y = centroid_coordinates[coord_index+1];
900
901      // Store x- and y- differentials for the edges of
902      // triangle k relative to the centroid
903      dxv0 = xv0 - x;
904      dxv1 = xv1 - x;
905      dxv2 = xv2 - x;
906      dyv0 = yv0 - y;
907      dyv1 = yv1 - y;
908      dyv2 = yv2 - y;
909      // Compute the minimum distance from the centroid to an edge
910      //demin=min(dxv0*dxv0 +dyv0*dyv0, min(dxv1*dxv1+dyv1*dyv1, dxv2*dxv2+dyv2*dyv2));
911      //demin=sqrt(demin);
912    }
913
914
915
916    if (number_of_boundaries[k]<=1)
917    {
918      //==============================================
919      // Number of boundaries <= 1
920      //==============================================
921
922
923      // If no boundaries, auxiliary triangle is formed
924      // from the centroids of the three neighbours
925      // If one boundary, auxiliary triangle is formed
926      // from this centroid and its two neighbours
927
928      k0 = surrogate_neighbours[k3];
929      k1 = surrogate_neighbours[k3 + 1];
930      k2 = surrogate_neighbours[k3 + 2];
931
932      // Test to see whether we accept the surrogate neighbours
933      // Note that if ki is replaced with k in more than 1 neighbour, then the
934      // triangle area will be zero, and a first order extrapolation will be
935      // used
936      if(stage_centroid_values[k2]<=elevation_centroid_values[k2]){
937          k2 = k ;
938      }
939      if(stage_centroid_values[k0]<=elevation_centroid_values[k0]){
940          k0 = k ;
941      }
942      if(stage_centroid_values[k1]<=elevation_centroid_values[k1]){
943          k1 = k ;
944      }
945      // Alternative approach (BRUTAL) -- if the max neighbour bed elevation is greater
946      // than the min neighbour stage, then we use first order extrapolation
947      //bedmax = max(elevation_centroid_values[k],
948      //             max(elevation_centroid_values[k0],
949      //                 max(elevation_centroid_values[k1], elevation_centroid_values[k2])));
950      //stagemin = min(stage_centroid_values[k],
951      //             min(stage_centroid_values[k0],
952      //                 min(stage_centroid_values[k1], stage_centroid_values[k2])));
953      //
954      //if(stagemin < bedmax){
955      //   // This will cause first order extrapolation
956      //   k2 = k;
957      //   k0 = k;
958      //   k1 = k;
959      //}
960
961      // Get the auxiliary triangle's vertex coordinates
962      // (really the centroids of neighbouring triangles)
963      coord_index = 2*k0;
964      x0 = centroid_coordinates[coord_index];
965      y0 = centroid_coordinates[coord_index+1];
966
967      coord_index = 2*k1;
968      x1 = centroid_coordinates[coord_index];
969      y1 = centroid_coordinates[coord_index+1];
970
971      coord_index = 2*k2;
972      x2 = centroid_coordinates[coord_index];
973      y2 = centroid_coordinates[coord_index+1];
974
975      // compute the maximum distance from the centroid to a neighbouring
976      // centroid
977      //dcmax=max( (x0-x)*(x0-x) + (y0-y)*(y0-y),
978      //           max((x1-x)*(x1-x) + (y1-y)*(y1-y),
979      //               (x2-x)*(x2-x) + (y2-y)*(y2-y)));
980      //dcmax=sqrt(dcmax);
981      //// Ratio of centroid to edge distance -- useful in attempting to adapt limiter
982      //if(dcmax>0.){
983      //    r0scale=demin/dcmax;
984      //    //printf("%f \n", r0scale);
985      //}else{
986      //    r0scale=0.5;
987      //}
988
989      // Store x- and y- differentials for the vertices
990      // of the auxiliary triangle
991      dx1 = x1 - x0;
992      dx2 = x2 - x0;
993      dy1 = y1 - y0;
994      dy2 = y2 - y0;
995
996      // Calculate 2*area of the auxiliary triangle
997      // The triangle is guaranteed to be counter-clockwise
998      area2 = dy2*dx1 - dy1*dx2;
999
1000      // If the mesh is 'weird' near the boundary,
1001      // the triangle might be flat or clockwise
1002      // Default to zero gradient
1003      if (area2 <= 0)
1004      {
1005          //printf("Error negative triangle area \n");
1006          //report_python_error(AT, "Negative triangle area");
1007          //return -1;
1008
1009          stage_edge_values[k3]   = stage_centroid_values[k];
1010          stage_edge_values[k3+1] = stage_centroid_values[k];
1011          stage_edge_values[k3+2] = stage_centroid_values[k];
1012          xmom_edge_values[k3]    = xmom_centroid_values[k];
1013          xmom_edge_values[k3+1]  = xmom_centroid_values[k];
1014          xmom_edge_values[k3+2]  = xmom_centroid_values[k];
1015          ymom_edge_values[k3]    = ymom_centroid_values[k];
1016          ymom_edge_values[k3+1]  = ymom_centroid_values[k];
1017          ymom_edge_values[k3+2]  = ymom_centroid_values[k];
1018
1019          continue;
1020      }
1021
1022      // Calculate heights of neighbouring cells
1023      hc = stage_centroid_values[k]  - elevation_centroid_values[k];
1024      h0 = stage_centroid_values[k0] - elevation_centroid_values[k0];
1025      h1 = stage_centroid_values[k1] - elevation_centroid_values[k1];
1026      h2 = stage_centroid_values[k2] - elevation_centroid_values[k2];
1027      hmin = min(min(h0, min(h1, h2)), hc);
1028      //hmin = min(h0, min(h1, h2));
1029      //hmin = max(hmin, 0.0);
1030      //hfactor = hc/(hc + 1.0);
1031
1032      hfactor = 0.0;
1033      //if (hmin > 0.001)
1034      if (hmin > 0.)
1035      //if (hc>0.0)
1036      {
1037        hfactor = 1.0 ;//hmin/(hmin + 0.004);
1038        //hfactor=hmin/(hmin + 0.004);
1039      }
1040
1041      if (optimise_dry_cells)
1042      {
1043        // Check if linear reconstruction is necessary for triangle k
1044        // This check will exclude dry cells.
1045
1046        //hmax = max(h0, max(h1, max(h2, hc)));
1047        hmax = max(h0, max(h1, h2));
1048        if (hmax < epsilon)
1049        {
1050            continue;
1051        }
1052      }
1053
1054      //-----------------------------------
1055      // stage
1056      //-----------------------------------
1057
1058      // Calculate the difference between vertex 0 of the auxiliary
1059      // triangle and the centroid of triangle k
1060      dq0 = stage_centroid_values[k0] - stage_centroid_values[k];
1061
1062      // Calculate differentials between the vertices
1063      // of the auxiliary triangle (centroids of neighbouring triangles)
1064      dq1 = stage_centroid_values[k1] - stage_centroid_values[k0];
1065      dq2 = stage_centroid_values[k2] - stage_centroid_values[k0];
1066
1067      inv_area2 = 1.0/area2;
1068      // Calculate the gradient of stage on the auxiliary triangle
1069      a = dy2*dq1 - dy1*dq2;
1070      a *= inv_area2;
1071      b = dx1*dq2 - dx2*dq1;
1072      b *= inv_area2;
1073
1074      // Calculate provisional jumps in stage from the centroid
1075      // of triangle k to its vertices, to be limited
1076      dqv[0] = a*dxv0 + b*dyv0;
1077      dqv[1] = a*dxv1 + b*dyv1;
1078      dqv[2] = a*dxv2 + b*dyv2;
1079
1080      // Now we want to find min and max of the centroid and the
1081      // vertices of the auxiliary triangle and compute jumps
1082      // from the centroid to the min and max
1083      find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
1084
1085      beta_tmp = beta_w_dry + (beta_w - beta_w_dry) * hfactor;
1086
1090
1091      //for (i=0;i<3;i++)
1092      stage_edge_values[k3+0] = stage_centroid_values[k] + dqv[0];
1093      stage_edge_values[k3+1] = stage_centroid_values[k] + dqv[1];
1094      stage_edge_values[k3+2] = stage_centroid_values[k] + dqv[2];
1095
1096      //-----------------------------------
1097      // xmomentum
1098      //-----------------------------------
1099
1100      // Calculate the difference between vertex 0 of the auxiliary
1101      // triangle and the centroid of triangle k
1102      dq0 = xmom_centroid_values[k0] - xmom_centroid_values[k];
1103
1104      // Calculate differentials between the vertices
1105      // of the auxiliary triangle
1106      dq1 = xmom_centroid_values[k1] - xmom_centroid_values[k0];
1107      dq2 = xmom_centroid_values[k2] - xmom_centroid_values[k0];
1108
1109      // Calculate the gradient of xmom on the auxiliary triangle
1110      a = dy2*dq1 - dy1*dq2;
1111      a *= inv_area2;
1112      b = dx1*dq2 - dx2*dq1;
1113      b *= inv_area2;
1114
1115      // Calculate provisional jumps in stage from the centroid
1116      // of triangle k to its vertices, to be limited
1117      dqv[0] = a*dxv0+b*dyv0;
1118      dqv[1] = a*dxv1+b*dyv1;
1119      dqv[2] = a*dxv2+b*dyv2;
1120
1121      // Now we want to find min and max of the centroid and the
1122      // vertices of the auxiliary triangle and compute jumps
1123      // from the centroid to the min and max
1124      //
1125      find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
1126
1127      beta_tmp = beta_uh_dry + (beta_uh - beta_uh_dry) * hfactor;
1128
1132
1133
1134      for (i=0; i < 3; i++)
1135      {
1136        xmom_edge_values[k3+i] = xmom_centroid_values[k] + dqv[i];
1137      }
1138
1139      //-----------------------------------
1140      // ymomentum
1141      //-----------------------------------
1142
1143      // Calculate the difference between vertex 0 of the auxiliary
1144      // triangle and the centroid of triangle k
1145      dq0 = ymom_centroid_values[k0] - ymom_centroid_values[k];
1146
1147      // Calculate differentials between the vertices
1148      // of the auxiliary triangle
1149      dq1 = ymom_centroid_values[k1] - ymom_centroid_values[k0];
1150      dq2 = ymom_centroid_values[k2] - ymom_centroid_values[k0];
1151
1152      // Calculate the gradient of xmom on the auxiliary triangle
1153      a = dy2*dq1 - dy1*dq2;
1154      a *= inv_area2;
1155      b = dx1*dq2 - dx2*dq1;
1156      b *= inv_area2;
1157
1158      // Calculate provisional jumps in stage from the centroid
1159      // of triangle k to its vertices, to be limited
1160      dqv[0] = a*dxv0 + b*dyv0;
1161      dqv[1] = a*dxv1 + b*dyv1;
1162      dqv[2] = a*dxv2 + b*dyv2;
1163
1164      // Now we want to find min and max of the centroid and the
1165      // vertices of the auxiliary triangle and compute jumps
1166      // from the centroid to the min and max
1167      //
1168      find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
1169
1170      beta_tmp = beta_vh_dry + (beta_vh - beta_vh_dry) * hfactor;
1171
1175
1176      for (i=0;i<3;i++)
1177      {
1178        ymom_edge_values[k3 + i] = ymom_centroid_values[k] + dqv[i];
1179      }
1180
1181    } // End number_of_boundaries <=1
1182    else
1183    {
1184
1185      //==============================================
1186      // Number of boundaries == 2
1187      //==============================================
1188
1189      // One internal neighbour and gradient is in direction of the neighbour's centroid
1190
1191      // Find the only internal neighbour (k1?)
1192      for (k2 = k3; k2 < k3 + 3; k2++)
1193      {
1194      // Find internal neighbour of triangle k
1195      // k2 indexes the edges of triangle k
1196
1197          if (surrogate_neighbours[k2] != k)
1198          {
1199             break;
1200          }
1201      }
1202
1203      if ((k2 == k3 + 3))
1204      {
1205        // If we didn't find an internal neighbour
1207        return -1;
1208      }
1209
1210      k1 = surrogate_neighbours[k2];
1211
1212      // The coordinates of the triangle are already (x,y).
1213      // Get centroid of the neighbour (x1,y1)
1214      coord_index = 2*k1;
1215      x1 = centroid_coordinates[coord_index];
1216      y1 = centroid_coordinates[coord_index + 1];
1217
1218      // Compute x- and y- distances between the centroid of
1219      // triangle k and that of its neighbour
1220      dx1 = x1 - x;
1221      dy1 = y1 - y;
1222
1223      // Set area2 as the square of the distance
1224      area2 = dx1*dx1 + dy1*dy1;
1225
1226      // Set dx2=(x1-x0)/((x1-x0)^2+(y1-y0)^2)
1227      // and dy2=(y1-y0)/((x1-x0)^2+(y1-y0)^2) which
1228      // respectively correspond to the x- and y- gradients
1229      // of the conserved quantities
1230      dx2 = 1.0/area2;
1231      dy2 = dx2*dy1;
1232      dx2 *= dx1;
1233
1234
1235      //-----------------------------------
1236      // stage
1237      //-----------------------------------
1238
1239      // Compute differentials
1240      dq1 = stage_centroid_values[k1] - stage_centroid_values[k];
1241
1242      // Calculate the gradient between the centroid of triangle k
1243      // and that of its neighbour
1244      a = dq1*dx2;
1245      b = dq1*dy2;
1246
1247      // Calculate provisional edge jumps, to be limited
1248      dqv[0] = a*dxv0 + b*dyv0;
1249      dqv[1] = a*dxv1 + b*dyv1;
1250      dqv[2] = a*dxv2 + b*dyv2;
1251
1252      // Now limit the jumps
1253      if (dq1>=0.0)
1254      {
1255        qmin=0.0;
1256        qmax=dq1;
1257      }
1258      else
1259      {
1260        qmin = dq1;
1261        qmax = 0.0;
1262      }
1263
1266
1267      //for (i=0; i < 3; i++)
1268      //{
1269      stage_edge_values[k3] = stage_centroid_values[k] + dqv[0];
1270      stage_edge_values[k3 + 1] = stage_centroid_values[k] + dqv[1];
1271      stage_edge_values[k3 + 2] = stage_centroid_values[k] + dqv[2];
1272      //}
1273
1274      //-----------------------------------
1275      // xmomentum
1276      //-----------------------------------
1277
1278      // Compute differentials
1279      dq1 = xmom_centroid_values[k1] - xmom_centroid_values[k];
1280
1281      // Calculate the gradient between the centroid of triangle k
1282      // and that of its neighbour
1283      a = dq1*dx2;
1284      b = dq1*dy2;
1285
1286      // Calculate provisional edge jumps, to be limited
1287      dqv[0] = a*dxv0+b*dyv0;
1288      dqv[1] = a*dxv1+b*dyv1;
1289      dqv[2] = a*dxv2+b*dyv2;
1290
1291      // Now limit the jumps
1292      if (dq1 >= 0.0)
1293      {
1294        qmin = 0.0;
1295        qmax = dq1;
1296      }
1297      else
1298      {
1299        qmin = dq1;
1300        qmax = 0.0;
1301      }
1302
1305
1306      //for (i=0;i<3;i++)
1307      //xmom_edge_values[k3] = xmom_centroid_values[k] + dqv[0];
1308      //xmom_edge_values[k3 + 1] = xmom_centroid_values[k] + dqv[1];
1309      //xmom_edge_values[k3 + 2] = xmom_centroid_values[k] + dqv[2];
1310
1311      for (i = 0; i < 3;i++)
1312      {
1313          xmom_edge_values[k3 + i] = xmom_centroid_values[k] + dqv[i];
1314      }
1315
1316      //-----------------------------------
1317      // ymomentum
1318      //-----------------------------------
1319
1320      // Compute differentials
1321      dq1 = ymom_centroid_values[k1] - ymom_centroid_values[k];
1322
1323      // Calculate the gradient between the centroid of triangle k
1324      // and that of its neighbour
1325      a = dq1*dx2;
1326      b = dq1*dy2;
1327
1328      // Calculate provisional edge jumps, to be limited
1329      dqv[0] = a*dxv0 + b*dyv0;
1330      dqv[1] = a*dxv1 + b*dyv1;
1331      dqv[2] = a*dxv2 + b*dyv2;
1332
1333      // Now limit the jumps
1334      if (dq1>=0.0)
1335      {
1336        qmin = 0.0;
1337        qmax = dq1;
1338      }
1339      else
1340      {
1341        qmin = dq1;
1342        qmax = 0.0;
1343      }
1344
1347
1348      for (i=0;i<3;i++)
1349              {
1350              ymom_edge_values[k3 + i] = ymom_centroid_values[k] + dqv[i];
1351              }
1352    } // else [number_of_boundaries==2]
1353  } // for k=0 to number_of_elements-1
1354
1355  // Compute vertex values of quantities
1356  for (k=0; k<number_of_elements; k++){
1357      k3=3*k;
1358
1359      // Compute stage vertex values
1360      stage_vertex_values[k3] = stage_edge_values[k3+1] + stage_edge_values[k3+2] -stage_edge_values[k3] ;
1361      stage_vertex_values[k3+1] =  stage_edge_values[k3] + stage_edge_values[k3+2]-stage_edge_values[k3+1];
1362      stage_vertex_values[k3+2] =  stage_edge_values[k3] + stage_edge_values[k3+1]-stage_edge_values[k3+2];
1363
1364     // Compute xmom vertex values
1365      xmom_vertex_values[k3] = xmom_edge_values[k3+1] + xmom_edge_values[k3+2] -xmom_edge_values[k3] ;
1366      xmom_vertex_values[k3+1] =  xmom_edge_values[k3] + xmom_edge_values[k3+2]-xmom_edge_values[k3+1];
1367      xmom_vertex_values[k3+2] =  xmom_edge_values[k3] + xmom_edge_values[k3+1]-xmom_edge_values[k3+2];
1368
1369      // Compute ymom vertex values
1370      ymom_vertex_values[k3] = ymom_edge_values[k3+1] + ymom_edge_values[k3+2] -ymom_edge_values[k3] ;
1371      ymom_vertex_values[k3+1] =  ymom_edge_values[k3] + ymom_edge_values[k3+2]-ymom_edge_values[k3+1];
1372      ymom_vertex_values[k3+2] =  ymom_edge_values[k3] + ymom_edge_values[k3+1]-ymom_edge_values[k3+2];
1373
1374      // If needed, convert from velocity to momenta
1375      if(extrapolate_velocity_second_order==1){
1376          //Convert velocity back to momenta at centroids
1377          xmom_centroid_values[k] = xmom_centroid_store[k];
1378          ymom_centroid_values[k] = ymom_centroid_store[k];
1379
1380          // Re-compute momenta at edges
1381          for (i=0; i<3; i++){
1382              de[i] = max(stage_edge_values[k3+i]-elevation_edge_values[k3+i],0.0);
1383              xmom_edge_values[k3+i]=xmom_edge_values[k3+i]*de[i];
1384              ymom_edge_values[k3+i]=ymom_edge_values[k3+i]*de[i];
1385          }
1386
1387          // Re-compute momenta at vertices
1388          for (i=0; i<3; i++){
1389              de[i] = max(stage_vertex_values[k3+i]-elevation_vertex_values[k3+i],0.0);
1390              xmom_vertex_values[k3+i]=xmom_vertex_values[k3+i]*de[i];
1391              ymom_vertex_values[k3+i]=ymom_vertex_values[k3+i]*de[i];
1392          }
1393
1394      }
1395
1396
1397  }
1398
1399  free(xmom_centroid_store);
1400  free(ymom_centroid_store);
1401  free(stage_centroid_store);
1402
1403  return 0;
1404}
1405
1406//=========================================================================
1407// Python Glue
1408//=========================================================================
1409
1410
1411//========================================================================
1412// Compute fluxes
1413//========================================================================
1414
1415// Modified central flux function
1416
1417PyObject *compute_fluxes_ext_central(PyObject *self, PyObject *args) {
1418  /*Compute all fluxes and the timestep suitable for all volumes
1419    in domain.
1420
1421    Compute total flux for each conserved quantity using "flux_function_central"
1422
1423    Fluxes across each edge are scaled by edgelengths and summed up
1424    Resulting flux is then scaled by area and stored in
1425    explicit_update for each of the three conserved quantities
1426    stage, xmomentum and ymomentum
1427
1428    The maximal allowable speed computed by the flux_function for each volume
1429    is converted to a timestep that must not be exceeded. The minimum of
1430    those is computed as the next overall timestep.
1431
1432    Python call:
1433    domain.timestep = compute_fluxes(timestep,
1434                                     domain.epsilon,
1435                     domain.H0,
1436                                     domain.g,
1437                                     domain.neighbours,
1438                                     domain.neighbour_edges,
1439                                     domain.normals,
1440                                     domain.edgelengths,
1442                                     domain.areas,
1443                                     tri_full_flag,
1444                                     Stage.edge_values,
1445                                     Xmom.edge_values,
1446                                     Ymom.edge_values,
1447                                     Bed.edge_values,
1448                                     Stage.boundary_values,
1449                                     Xmom.boundary_values,
1450                                     Ymom.boundary_values,
1451                                     Stage.explicit_update,
1452                                     Xmom.explicit_update,
1453                                     Ymom.explicit_update,
1455                                     optimise_dry_cells,
1456                                     stage.centroid_values,
1457                                     bed.centroid_values)
1458
1459
1460    Post conditions:
1461      domain.explicit_update is reset to computed flux values
1462      domain.timestep is set to the largest step satisfying all volumes.
1463
1464
1465  */
1466
1467
1468  PyArrayObject *neighbours, *neighbour_edges,
1470    *tri_full_flag,
1471    *stage_edge_values,
1472    *xmom_edge_values,
1473    *ymom_edge_values,
1474    *bed_edge_values,
1475    *stage_boundary_values,
1476    *xmom_boundary_values,
1477    *ymom_boundary_values,
1478    *boundary_flux_type,
1479    *stage_explicit_update,
1480    *xmom_explicit_update,
1481    *ymom_explicit_update,
1482    *already_computed_flux, //Tracks whether the flux across an edge has already been computed
1483    *max_speed_array, //Keeps track of max speeds for each triangle
1484    *stage_centroid_values,
1485    *bed_centroid_values,
1486    *bed_vertex_values;
1487
1488  double timestep, epsilon, H0, g;
1489  int optimise_dry_cells;
1490
1491  // Convert Python arguments to C
1492  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOOiOOO",
1493            &timestep,
1494            &epsilon,
1495            &H0,
1496            &g,
1497            &neighbours,
1498            &neighbour_edges,
1499            &normals,
1501            &tri_full_flag,
1502            &stage_edge_values,
1503            &xmom_edge_values,
1504            &ymom_edge_values,
1505            &bed_edge_values,
1506            &stage_boundary_values,
1507            &xmom_boundary_values,
1508            &ymom_boundary_values,
1509            &boundary_flux_type,
1510            &stage_explicit_update,
1511            &xmom_explicit_update,
1512            &ymom_explicit_update,
1514            &max_speed_array,
1515            &optimise_dry_cells,
1516            &stage_centroid_values,
1517            &bed_centroid_values,
1518            &bed_vertex_values)) {
1519    report_python_error(AT, "could not parse input arguments");
1520    return NULL;
1521  }
1522
1523  // check that numpy array objects arrays are C contiguous memory
1524  CHECK_C_CONTIG(neighbours);
1525  CHECK_C_CONTIG(neighbour_edges);
1526  CHECK_C_CONTIG(normals);
1527  CHECK_C_CONTIG(edgelengths);
1529  CHECK_C_CONTIG(areas);
1530  CHECK_C_CONTIG(tri_full_flag);
1531  CHECK_C_CONTIG(stage_edge_values);
1532  CHECK_C_CONTIG(xmom_edge_values);
1533  CHECK_C_CONTIG(ymom_edge_values);
1534  CHECK_C_CONTIG(bed_edge_values);
1535  CHECK_C_CONTIG(stage_boundary_values);
1536  CHECK_C_CONTIG(xmom_boundary_values);
1537  CHECK_C_CONTIG(ymom_boundary_values);
1538  CHECK_C_CONTIG(boundary_flux_type);
1539  CHECK_C_CONTIG(stage_explicit_update);
1540  CHECK_C_CONTIG(xmom_explicit_update);
1541  CHECK_C_CONTIG(ymom_explicit_update);
1543  CHECK_C_CONTIG(max_speed_array);
1544  CHECK_C_CONTIG(stage_centroid_values);
1545  CHECK_C_CONTIG(bed_centroid_values);
1546  CHECK_C_CONTIG(bed_vertex_values);
1547
1548  int number_of_elements = stage_edge_values -> dimensions[0];
1549
1550  // Call underlying flux computation routine and update
1551  // the explicit update arrays
1552  timestep = _compute_fluxes_central(number_of_elements,
1553                     timestep,
1554                     epsilon,
1555                     H0,
1556                     g,
1557                     (long*) neighbours -> data,
1558                     (long*) neighbour_edges -> data,
1559                     (double*) normals -> data,
1560                     (double*) edgelengths -> data,
1562                     (double*) areas -> data,
1563                     (long*) tri_full_flag -> data,
1564                     (double*) stage_edge_values -> data,
1565                     (double*) xmom_edge_values -> data,
1566                     (double*) ymom_edge_values -> data,
1567                     (double*) bed_edge_values -> data,
1568                     (double*) stage_boundary_values -> data,
1569                     (double*) xmom_boundary_values -> data,
1570                     (double*) ymom_boundary_values -> data,
1571                     (long*)   boundary_flux_type -> data,
1572                     (double*) stage_explicit_update -> data,
1573                     (double*) xmom_explicit_update -> data,
1574                     (double*) ymom_explicit_update -> data,
1576                     (double*) max_speed_array -> data,
1577                     optimise_dry_cells,
1578                     (double*) stage_centroid_values -> data,
1579                     (double*) bed_centroid_values -> data,
1580                     (double*) bed_vertex_values -> data);
1581
1582  // Return updated flux timestep
1583  return Py_BuildValue("d", timestep);
1584}
1585
1586
1587PyObject *flux_function_central(PyObject *self, PyObject *args) {
1588  //
1589  // Gateway to innermost flux function.
1590  // This is only used by the unit tests as the c implementation is
1591  // normally called by compute_fluxes in this module.
1592
1593
1594  PyArrayObject *normal, *ql, *qr,  *edgeflux;
1595  double g, epsilon, max_speed, H0, zl, zr;
1596  double h0, limiting_threshold, pressure_flux;
1597
1598  if (!PyArg_ParseTuple(args, "OOOddOddd",
1599            &normal, &ql, &qr, &zl, &zr, &edgeflux,
1600            &epsilon, &g, &H0)) {
1601      report_python_error(AT, "could not parse input arguments");
1602      return NULL;
1603  }
1604
1605
1606  h0 = H0*H0; // This ensures a good balance when h approaches H0.
1607              // But evidence suggests that h0 can be as little as
1608              // epsilon!
1609
1610  limiting_threshold = 10*H0; // Avoid applying limiter below this
1611                              // threshold for performance reasons.
1612                              // See ANUGA manual under flux limiting
1613
1614  pressure_flux = 0.0; // Store the water pressure related component of the flux
1615  _flux_function_central((double*) ql -> data,
1616                         (double*) qr -> data,
1617                         zl,
1618                         zr,
1619                         ((double*) normal -> data)[0],
1620                         ((double*) normal -> data)[1],
1621                         epsilon, h0, limiting_threshold,
1622                         g,
1623                         (double*) edgeflux -> data,
1624                         &max_speed,
1625             &pressure_flux);
1626
1627  return Py_BuildValue("d", max_speed);
1628}
1629
1630//========================================================================
1631// Gravity
1632//========================================================================
1633
1634PyObject *gravity(PyObject *self, PyObject *args) {
1635  //
1636  //  gravity(g, h, v, x, xmom, ymom)
1637  //
1638
1639
1640  PyArrayObject *h, *v, *x, *xmom, *ymom;
1641  int k, N, k3, k6;
1642  double g, avg_h, zx, zy;
1643  double x0, y0, x1, y1, x2, y2, z0, z1, z2;
1644  //double epsilon;
1645
1646  if (!PyArg_ParseTuple(args, "dOOOOO",
1647                        &g, &h, &v, &x,
1648                        &xmom, &ymom)) {
1649    //&epsilon)) {
1650    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: gravity could not parse input arguments");
1651    return NULL;
1652  }
1653
1654  // check that numpy array objects arrays are C contiguous memory
1655  CHECK_C_CONTIG(h);
1656  CHECK_C_CONTIG(v);
1657  CHECK_C_CONTIG(x);
1658  CHECK_C_CONTIG(xmom);
1659  CHECK_C_CONTIG(ymom);
1660
1661  N = h -> dimensions[0];
1662  for (k=0; k<N; k++) {
1663    k3 = 3*k;  // base index
1664
1665    // Get bathymetry
1666    z0 = ((double*) v -> data)[k3 + 0];
1667    z1 = ((double*) v -> data)[k3 + 1];
1668    z2 = ((double*) v -> data)[k3 + 2];
1669
1670    // Optimise for flat bed
1671    // Note (Ole): This didn't produce measurable speed up.
1672    // Revisit later
1673    //if (fabs(z0-z1)<epsilon && fabs(z1-z2)<epsilon) {
1674    //  continue;
1675    //}
1676
1677    // Get average depth from centroid values
1678    avg_h = ((double *) h -> data)[k];
1679
1680    // Compute bed slope
1681    k6 = 6*k;  // base index
1682
1683    x0 = ((double*) x -> data)[k6 + 0];
1684    y0 = ((double*) x -> data)[k6 + 1];
1685    x1 = ((double*) x -> data)[k6 + 2];
1686    y1 = ((double*) x -> data)[k6 + 3];
1687    x2 = ((double*) x -> data)[k6 + 4];
1688    y2 = ((double*) x -> data)[k6 + 5];
1689
1690
1691    _gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2, &zx, &zy);
1692
1693    // Update momentum
1694    ((double*) xmom -> data)[k] += -g*zx*avg_h;
1695    ((double*) ymom -> data)[k] += -g*zy*avg_h;
1696  }
1697
1698  return Py_BuildValue("");
1699}
1700
1701
1702PyObject *extrapolate_second_order_edge_sw(PyObject *self, PyObject *args) {
1703  /*Compute the edge values based on a linear reconstruction
1704    on each triangle
1705
1706    Python call:
1707    extrapolate_second_order_sw(domain.surrogate_neighbours,
1708                                domain.number_of_boundaries
1709                                domain.centroid_coordinates,
1710                                Stage.centroid_values
1711                                Xmom.centroid_values
1712                                Ymom.centroid_values
1713                                domain.edge_coordinates,
1714                                Stage.edge_values,
1715                                Xmom.edge_values,
1716                                Ymom.edge_values)
1717
1718    Post conditions:
1719        The edges of each triangle have values from a
1720        limited linear reconstruction
1721        based on centroid values
1722
1723  */
1724  PyArrayObject *surrogate_neighbours,
1725    *number_of_boundaries,
1726    *centroid_coordinates,
1727    *stage_centroid_values,
1728    *xmom_centroid_values,
1729    *ymom_centroid_values,
1730    *elevation_centroid_values,
1731    *edge_coordinates,
1732    *stage_edge_values,
1733    *xmom_edge_values,
1734    *ymom_edge_values,
1735    *elevation_edge_values,
1736    *stage_vertex_values,
1737    *xmom_vertex_values,
1738    *ymom_vertex_values,
1739    *elevation_vertex_values;
1740
1741  PyObject *domain;
1742
1743
1744  double beta_w, beta_w_dry, beta_uh, beta_uh_dry, beta_vh, beta_vh_dry;
1745  double minimum_allowed_height, epsilon;
1746  int optimise_dry_cells, number_of_elements, extrapolate_velocity_second_order, e, e2;
1747
1748  // Convert Python arguments to C
1749  if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOOOOOii",
1750            &domain,
1751            &surrogate_neighbours,
1752            &number_of_boundaries,
1753            &centroid_coordinates,
1754            &stage_centroid_values,
1755            &xmom_centroid_values,
1756            &ymom_centroid_values,
1757            &elevation_centroid_values,
1758            &edge_coordinates,
1759            &stage_edge_values,
1760            &xmom_edge_values,
1761            &ymom_edge_values,
1762            &elevation_edge_values,
1763            &stage_vertex_values,
1764            &xmom_vertex_values,
1765            &ymom_vertex_values,
1766            &elevation_vertex_values,
1767            &optimise_dry_cells,
1768            &extrapolate_velocity_second_order)) {
1769
1770    report_python_error(AT, "could not parse input arguments");
1771    return NULL;
1772  }
1773
1774  // check that numpy array objects arrays are C contiguous memory
1775  CHECK_C_CONTIG(surrogate_neighbours);
1776  CHECK_C_CONTIG(number_of_boundaries);
1777  CHECK_C_CONTIG(centroid_coordinates);
1778  CHECK_C_CONTIG(stage_centroid_values);
1779  CHECK_C_CONTIG(xmom_centroid_values);
1780  CHECK_C_CONTIG(ymom_centroid_values);
1781  CHECK_C_CONTIG(elevation_centroid_values);
1782  CHECK_C_CONTIG(edge_coordinates);
1783  CHECK_C_CONTIG(stage_edge_values);
1784  CHECK_C_CONTIG(xmom_edge_values);
1785  CHECK_C_CONTIG(ymom_edge_values);
1786  CHECK_C_CONTIG(elevation_edge_values);
1787  CHECK_C_CONTIG(stage_vertex_values);
1788  CHECK_C_CONTIG(xmom_vertex_values);
1789  CHECK_C_CONTIG(ymom_vertex_values);
1790  CHECK_C_CONTIG(elevation_vertex_values);
1791
1792  // Get the safety factor beta_w, set in the config.py file.
1793  // This is used in the limiting process
1794
1795
1796  beta_w                 = get_python_double(domain,"beta_w");
1797  beta_w_dry             = get_python_double(domain,"beta_w_dry");
1798  beta_uh                = get_python_double(domain,"beta_uh");
1799  beta_uh_dry            = get_python_double(domain,"beta_uh_dry");
1800  beta_vh                = get_python_double(domain,"beta_vh");
1801  beta_vh_dry            = get_python_double(domain,"beta_vh_dry");
1802
1803  minimum_allowed_height = get_python_double(domain,"minimum_allowed_height");
1804  epsilon                = get_python_double(domain,"epsilon");
1805
1806  number_of_elements = stage_centroid_values -> dimensions[0];
1807
1808  //printf("In C before Extrapolate");
1809  //e=1;
1810  // Call underlying computational routine
1811  e = _extrapolate_second_order_edge_sw(number_of_elements,
1812                   epsilon,
1813                   minimum_allowed_height,
1814                   beta_w,
1815                   beta_w_dry,
1816                   beta_uh,
1817                   beta_uh_dry,
1818                   beta_vh,
1819                   beta_vh_dry,
1820                   (long*) surrogate_neighbours -> data,
1821                   (long*) number_of_boundaries -> data,
1822                   (double*) centroid_coordinates -> data,
1823                   (double*) stage_centroid_values -> data,
1824                   (double*) xmom_centroid_values -> data,
1825                   (double*) ymom_centroid_values -> data,
1826                   (double*) elevation_centroid_values -> data,
1827                   (double*) edge_coordinates -> data,
1828                   (double*) stage_edge_values -> data,
1829                   (double*) xmom_edge_values -> data,
1830                   (double*) ymom_edge_values -> data,
1831                   (double*) elevation_edge_values -> data,
1832                   (double*) stage_vertex_values -> data,
1833                   (double*) xmom_vertex_values -> data,
1834                   (double*) ymom_vertex_values -> data,
1835                   (double*) elevation_vertex_values -> data,
1836                   optimise_dry_cells,
1837                   extrapolate_velocity_second_order);
1838
1839  if (e == -1) {
1840    // Use error string set inside computational routine
1841    return NULL;
1842  }
1843
1844
1845  return Py_BuildValue("");
1846
1847}// extrapolate_second-order_edge_sw
1848//========================================================================
1849// Protect -- to prevent the water level from falling below the minimum
1850// bed_edge_value
1851//========================================================================
1852
1853PyObject *protect(PyObject *self, PyObject *args) {
1854  //
1855  //    protect(minimum_allowed_height, maximum_allowed_speed, wc, zc, xmomc, ymomc)
1856
1857
1858  PyArrayObject
1859  *wc,            //Stage at centroids
1860  *wv,            //Stage at vertices
1861  *zc,            //Elevation at centroids
1862  *zv,            //Elevation at vertices
1863  *xmomc,         //Momentums at centroids
1864  *ymomc;
1865
1866
1867  int N;
1868  double minimum_allowed_height, maximum_allowed_speed, epsilon;
1869
1870  // Convert Python arguments to C
1871  if (!PyArg_ParseTuple(args, "dddOOOOOO",
1872            &minimum_allowed_height,
1873            &maximum_allowed_speed,
1874            &epsilon,
1875            &wc, &wv, &zc, &zv, &xmomc, &ymomc)) {
1876    report_python_error(AT, "could not parse input arguments");
1877    return NULL;
1878  }
1879
1880  N = wc -> dimensions[0];
1881
1882  _protect(N,
1883       minimum_allowed_height,
1884       maximum_allowed_speed,
1885       epsilon,
1886       (double*) wc -> data,
1887       (double*) wv -> data,
1888       (double*) zc -> data,
1889       (double*) zv -> data,
1890       (double*) xmomc -> data,
1891       (double*) ymomc -> data);
1892
1893  return Py_BuildValue("");
1894}
1895
1896//========================================================================
1897// Method table for python module
1898//========================================================================
1899
1900static struct PyMethodDef MethodTable[] = {
1901  /* The cast of the function is necessary since PyCFunction values
1902   * only take two PyObject* parameters, and rotate() takes
1903   * three.
1904   */
1905  //{"rotate", (PyCFunction)rotate, METH_VARARGS | METH_KEYWORDS, "Print out"},
1906  {"compute_fluxes_ext_central", compute_fluxes_ext_central, METH_VARARGS, "Print out"},
1907  {"gravity_c",        gravity,            METH_VARARGS, "Print out"},
1908  {"flux_function_central", flux_function_central, METH_VARARGS, "Print out"},
1909  {"extrapolate_second_order_edge_sw", extrapolate_second_order_edge_sw, METH_VARARGS, "Print out"},
1910  {"protect",          protect, METH_VARARGS | METH_KEYWORDS, "Print out"},
1911  {NULL, NULL, 0, NULL}
1912};
1913
1914// Module initialisation
1915void initswb2_domain_ext(void){
1916  Py_InitModule("swb2_domain_ext", MethodTable);
1917
1918  import_array(); // Necessary for handling of NumPY structures
1919}
Note: See TracBrowser for help on using the repository browser.