# source:trunk/anuga_work/development/gareth/balanced_basic/swb2_domain_ext.c@8302

Last change on this file since 8302 was 8302, checked in by davies, 12 years ago

File size: 39.4 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 = 0.0;  // Could have been negative
70      u = 0.0;
71    } else {
72      u = *uh/(*h + h0/ *h);
73    }
74
75
76    // Adjust momentum to be consistent with speed
77    *uh = u * *h;
78  } else {
79    // We are in deep water - no need for limiting
80    u = *uh/ *h;
81  }
82
83  return u;
84}
85
86// Innermost flux function (using stage w=z+h)
87int _flux_function_central(double *q_left, double *q_right,
88                           double z_left, double z_right,
89                           double n1, double n2,
90                           double epsilon,
91                           double h0,
92                           double limiting_threshold,
93                           double g,
94                           double *edgeflux, double *max_speed,
95                           double *pressure_flux)
96{
97
98  /*Compute fluxes between volumes for the shallow water wave equation
99    cast in terms of the 'stage', w = h+z using
100    the 'central scheme' as described in
101
102    Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
103    Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
104    Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.
105
106    The implemented formula is given in equation (3.15) on page 714
107  */
108
109  int i;
110
111  double w_left, h_left, uh_left, vh_left, u_left;
112  double w_right, h_right, uh_right, vh_right, u_right;
113  double s_min, s_max, soundspeed_left, soundspeed_right;
114  double denom, inverse_denominator, z;
115  double uint, t1, t2, t3;
116  // Workspace (allocate once, use many)
117  static double q_left_rotated[3], q_right_rotated[3], flux_right[3], flux_left[3];
118
119  // Copy conserved quantities to protect from modification
120  q_left_rotated[0] = q_left[0];
121  q_right_rotated[0] = q_right[0];
122  q_left_rotated[1] = q_left[1];
123  q_right_rotated[1] = q_right[1];
124  q_left_rotated[2] = q_left[2];
125  q_right_rotated[2] = q_right[2];
126
127  // Align x- and y-momentum with x-axis
128  _rotate(q_left_rotated, n1, n2);
129  _rotate(q_right_rotated, n1, n2);
130
131  z = 0.5*(z_left + z_right); // Average elevation values.
132                              // Even though this will nominally allow
133                              // for discontinuities in the elevation data,
134                              // there is currently no numerical support for
135                              // this so results may be strange near
136                              // jumps in the bed.
137
138  // Compute speeds in x-direction
139  w_left = q_left_rotated[0];
140  h_left = w_left - z;
141  uh_left = q_left_rotated[1];
142  u_left = _compute_speed(&uh_left, &h_left,
143                          epsilon, h0, limiting_threshold);
144
145  w_right = q_right_rotated[0];
146  h_right = w_right - z;
147  uh_right = q_right_rotated[1];
148  u_right = _compute_speed(&uh_right, &h_right,
149                           epsilon, h0, limiting_threshold);
150
151  // Momentum in y-direction
152  vh_left  = q_left_rotated[2];
153  vh_right = q_right_rotated[2];
154
155  // Limit y-momentum if necessary
156  // Leaving this out, improves speed significantly (Ole 27/5/2009)
157  // All validation tests pass, so do we really need it anymore?
158  _compute_speed(&vh_left, &h_left,
159                 epsilon, h0, limiting_threshold);
160  _compute_speed(&vh_right, &h_right,
161                 epsilon, h0, limiting_threshold);
162
163  // Maximal and minimal wave speeds
164  soundspeed_left  = sqrt(g*h_left);
165  soundspeed_right = sqrt(g*h_right);
166
167  // Code to use fast square root optimisation if desired.
168  // Timings on AMD 64 for the Okushiri profile gave the following timings
169  //
170  // SQRT           Total    Flux
171  //=============================
172  //
173  // Ref            405s     152s
174  // Fast (dbl)     453s     173s
175  // Fast (sng)     437s     171s
176  //
177  // Consequently, there is currently (14/5/2009) no reason to use this
178  // approximation.
179
180  //soundspeed_left  = fast_squareroot_approximation(g*h_left);
181  //soundspeed_right = fast_squareroot_approximation(g*h_right);
182
183  s_max = max(u_left + soundspeed_left, u_right + soundspeed_right);
184  if (s_max < 0.0)
185  {
186    s_max = 0.0;
187  }
188
189  s_min = min(u_left - soundspeed_left, u_right - soundspeed_right);
190  if (s_min > 0.0)
191  {
192    s_min = 0.0;
193  }
194
195  // Flux formulas
196  flux_left[0] = u_left*h_left;
197  flux_left[1] = u_left*uh_left + 0.5*g*h_left*h_left;
198  flux_left[2] = u_left*vh_left;
199
200  flux_right[0] = u_right*h_right;
201  flux_right[1] = u_right*uh_right + 0.5*g*h_right*h_right;
202  flux_right[2] = u_right*vh_right;
203
204  // Flux computation
205  denom = s_max - s_min;
206  if (denom < epsilon)
207  { // FIXME (Ole): Try using h0 here
208    memset(edgeflux, 0, 3*sizeof(double));
209    //for(i=0;i<3;i++){
210    //    edgeflux[i] = _minmod(flux_left[i], flux_right[i]);
211    //}
212    *max_speed = 0.0;
213  }
214  else
215  {
216    inverse_denominator = 1.0/denom;
217    for (i = 0; i < 3; i++)
218    {
219      // Adjustment to the scheme by Kurganov and Lin 2007 Communications in Computational
220      // Physics 2:141-163
221      //uint = (s_max*q_right_rotated[i] - s_min*q_left_rotated[i] - (flux_right[i] - flux_left[i]))*inverse_denominator;
222      //t1 = (q_right_rotated[i] - uint);
223      //t2 = (-q_left_rotated[i] + uint);
224      //t3 = _minmod(t1,t2);
225      t3 = 0.0;
226      edgeflux[i] = s_max*flux_left[i] - s_min*flux_right[i];
227      edgeflux[i] += s_max*s_min*(q_right_rotated[i] - q_left_rotated[i] - t3);
228      //t1 = s_max*s_min*(q_right_rotated[i] - q_left_rotated[i] - t3);
229      //if(fabs(edgeflux[i])>fabs(t1)){
230      //      edgeflux[i]+=t1;
231      //}else{
232      //      edgeflux[i]+=fabs(edgeflux[i])*copysign(1.0,t1) ;
233      //}
234
235      edgeflux[i] *= inverse_denominator;
236    }
237
238    *pressure_flux = 0.5*g*( s_max*h_left*h_left -s_min*h_right*h_right)*inverse_denominator;
239    //if(edgeflux[0]<0.0){
240    //    //edgeflux[2] = edgeflux[0]*0.5*(vh_right/h_right + vh_left/h_left);
241    //    edgeflux[2] = edgeflux[0]*vh_right/h_right;
242    //}else{
243    //    edgeflux[2] = edgeflux[0]*vh_left/h_left;
244    //}
245
246    // Maximal wavespeed
247    *max_speed = max(fabs(s_max), fabs(s_min));
248
249    // Rotate back
250    _rotate(edgeflux, n1, -n2);
251  }
252
253  return 0;
254}
255
256
257// Computational function for flux computation
258double _compute_fluxes_central(int number_of_elements,
259        double timestep,
260        double epsilon,
261        double H0,
262        double g,
263        long* neighbours,
264        long* neighbour_edges,
265        double* normals,
266        double* edgelengths,
268        double* areas,
269        long* tri_full_flag,
270        double* stage_edge_values,
271        double* xmom_edge_values,
272        double* ymom_edge_values,
273        double* bed_edge_values,
274        double* stage_boundary_values,
275        double* xmom_boundary_values,
276        double* ymom_boundary_values,
277        double* stage_explicit_update,
278        double* xmom_explicit_update,
279        double* ymom_explicit_update,
281        double* max_speed_array,
282        int optimise_dry_cells,
283        double* stage_centroid_values,
284        double* bed_centroid_values,
285        double* bed_vertex_values) {
286    // Local variables
287    double max_speed, length, inv_area, zl, zr;
288    double h0 = H0*H0; // This ensures a good balance when h approaches H0.
289
290    double limiting_threshold = 10 * H0; // Avoid applying limiter below this
291    // threshold for performance reasons.
292    // See ANUGA manual under flux limiting
293    int k, i, m, n,j;
294    int ki, nm = 0, ki2; // Index shorthands
295
296    // Workspace (making them static actually made function slightly slower (Ole))
297    double ql[3], qr[3], edgeflux[3]; // Work array for summing up fluxes
298    double stage_edges[3];//Work array
299    double bedslope_work;
300    int neighbours_wet[3];//Work array
301    int useint;
302    double usedouble, scale_factor_shallow, bedtop, bedbot, pressure_flux;
303    //double depthtemp, max_bed_edge_value;
304    static long call = 1; // Static local variable flagging already computed flux
305
306    // Start computation
307    call++; // Flag 'id' of flux calculation for this timestep
308
309    // Set explicit_update to zero for all conserved_quantities.
310    // This assumes compute_fluxes called before forcing terms
311    memset((char*) stage_explicit_update, 0, number_of_elements * sizeof (double));
312    memset((char*) xmom_explicit_update, 0, number_of_elements * sizeof (double));
313    memset((char*) ymom_explicit_update, 0, number_of_elements * sizeof (double));
314
315    // For all triangles
316    for (k = 0; k < number_of_elements; k++) {
317        // Loop through neighbours and compute edge flux for each
318        for (i = 0; i < 3; i++) {
319            ki = k * 3 + i; // Linear index to edge i of triangle k
320
321            if (already_computed_flux[ki] == call) {
322                // We've already computed the flux across this edge
323                continue;
324            }
325
326            // Get left hand side values from triangle k, edge i
327            ql[0] = stage_edge_values[ki];
328            ql[1] = xmom_edge_values[ki];
329            ql[2] = ymom_edge_values[ki];
330            zl = bed_edge_values[ki];
331
332            // Get right hand side values either from neighbouring triangle
333            // or from boundary array (Quantities at neighbour on nearest face).
334            n = neighbours[ki];
335            if (n < 0) {
336                // Neighbour is a boundary condition
337                m = -n - 1; // Convert negative flag to boundary index
338
339                qr[0] = stage_boundary_values[m];
340                qr[1] = xmom_boundary_values[m];
341                qr[2] = ymom_boundary_values[m];
342                zr = zl; // Extend bed elevation to boundary
343            }
344            else {
345                // Neighbour is a real triangle
346                m = neighbour_edges[ki];
347                nm = n * 3 + m; // Linear index (triangle n, edge m)
348
349                qr[0] = stage_edge_values[nm];
350                qr[1] = xmom_edge_values[nm];
351                qr[2] = ymom_edge_values[nm];
352                zr = bed_edge_values[nm];
353            }
354
355            // Now we have values for this edge - both from left and right side.
356
357            if (optimise_dry_cells) {
358                // Check if flux calculation is necessary across this edge
359                // This check will exclude dry cells.
360                // This will also optimise cases where zl != zr as
361                // long as both are dry
362
363                if ((ql[0] - zl) < epsilon &&
364                        (qr[0] - zr) < epsilon) {
365                    // Cell boundary is dry
366
367                    already_computed_flux[ki] = call; // #k Done
368                    if (n >= 0) {
369                        already_computed_flux[nm] = call; // #n Done
370                    }
371
372                    max_speed = 0.0;
373                    continue;
374                }
375            }
376
377
378            if (fabs(zl-zr)>1.0e-10) {
379                report_python_error(AT,"Discontinuous Elevation");
380                return 0.0;
381            }
382
383            // If both centroids are dry, then the cells should not exchange flux
384            // NOTE: This also means no bedslope term -- which is arguably
385            // appropriate, since the depth is zero at the cell centroid
386            if( (stage_centroid_values[k]-bed_centroid_values[k]<=0.0)){
387                if(n>0){
388                    if(stage_centroid_values[n]-bed_centroid_values[n]<=0.0){
389                        continue;
390                    }
391                }else{
392                    // Treat boundaries -- if the boundary stage < zl, there is
393                    // no flux
394                    if(qr[0]< zl){
395                        continue;
396                    }
397                }
398            }
399
400            //If one centroid is dry, then extrapolate its edge values from the neighbouring centroid,
401            // unless the local centroid value is smaller
402            if(n>=0){
403                if((stage_centroid_values[k]-bed_centroid_values[k]<= 0.0)){
404                          ql[0]=max(min(qr[0],stage_centroid_values[k]),zl);
405                }
406                if(stage_centroid_values[n]-bed_centroid_values[n]<= 0.0){
407                        qr[0]=max(min(ql[0],stage_centroid_values[n]),zr);
408                }
409            }else{
410                // Treat the boundary case
411                if((stage_centroid_values[k]-bed_centroid_values[k]<=0.0)){
412                          ql[0]=max(min(qr[0],stage_centroid_values[k]),zl);
413                }
414
415            }
416
417            // Outward pointing normal vector (domain.normals[k, 2*i:2*i+2])
418            ki2 = 2 * ki; //k*6 + i*2
419
420            // Edge flux computation (triangle k, edge i)
421            _flux_function_central(ql, qr, zl, zr,
422                    normals[ki2], normals[ki2 + 1],
423                    epsilon, h0, limiting_threshold, g,
424                    edgeflux, &max_speed, &pressure_flux);
425            //_flux_function_fraccaro(ql, qr, zl, zr,
426            //        normals[ki2], normals[ki2 + 1],
427            //        epsilon, h0, limiting_threshold, g,
428            //        edgeflux, &max_speed);
429
430            //Find the bed vertex values associated with the edge
431            //if(i==0){
432            //    bedtop = max(bed_vertex_values[3*k+1], bed_vertex_values[3*k+2]);
433            //    bedbot = min(bed_vertex_values[3*k+1], bed_vertex_values[3*k+2]);
434            //}else{
435            //    if(i==1){
436            //    bedtop = max(bed_vertex_values[3*k+0], bed_vertex_values[3*k+2]);
437            //    bedbot = min(bed_vertex_values[3*k+0], bed_vertex_values[3*k+2]);
438            //    }else{
439            //    // i==2
440            //    bedtop = max(bed_vertex_values[3*k+1], bed_vertex_values[3*k+0]);
441            //    bedbot = min(bed_vertex_values[3*k+1], bed_vertex_values[3*k+0]);
442            //    }
443            //
444            //}
445
446            //if(fabs(bedtop-bedbot)<1.0e-10){
447            //    scale_factor_shallow = 1.0;
448            //}else{
449            //    scale_factor_shallow = min( (bedtop - bedbot)/(min(stage_centroid_values[k], bedtop) - bedbot),10.0);
450            //}
451
452            // Multiply edgeflux by edgelength
453            length = edgelengths[ki];
454            //length = edgelengths[ki]*scale_factor_shallow;
455            edgeflux[0] *= length;
456            edgeflux[1] *= length;
457            edgeflux[2] *= length;
458
459            // To correctly implement reflective boundary conditions would be something like this?
460            //if(n<0){
461            //   _rotate(edgeflux,normals[ki2],normals[ki2+1]);
462            //   edgeflux[1] = 0.0;
463            //   _rotate(edgeflux,-normals[ki2+1],normals[ki2]);
464            //}
465
466            // Update triangle k with flux from edge i
467            stage_explicit_update[k] -= edgeflux[0];
468            xmom_explicit_update[k] -= edgeflux[1];
469            ymom_explicit_update[k] -= edgeflux[2];
470
471            // Calculate the bed slope term, using the 'divergence form' from
472            // Alessandro Valiani and Lorenzo Begnudelli, Divergence Form for Bed Slope Source Term in Shallow Water Equations
473            // Journal of Hydraulic Engineering, 2006, 132:652-665
474            if((stage_centroid_values[k]>bed_centroid_values[k]+0.0)){
475                bedslope_work = -0.5*g*max(stage_centroid_values[k]-zl,0.)*(stage_centroid_values[k]-zl)*length;
476                xmom_explicit_update[k] -= normals[ki2]*bedslope_work;
477                ymom_explicit_update[k] -= normals[ki2+1]*bedslope_work;
478            }else{
479                // Treat nearly dry cells
480                //bedslope_work =-0.5*g*max(ql[0]-zl,0.)*(ql[0]-zl)*length; //
481                bedslope_work = -pressure_flux*length;
482                xmom_explicit_update[k] -= normals[ki2]*bedslope_work;
483                ymom_explicit_update[k] -= normals[ki2+1]*bedslope_work;
484            }
485                // The triangle is (partially?) dry
486                // Find the average stage of the wet neighbours edges -- set this as the within triangle stage for the bedslope term
487                // FIXME: Note that we keep recomputing this -- could improve efficiency
488                //for(j=0;j<3;j++){
489                //    // For each edge, check if each neighbouring triangle is wet
490                //    useint = neighbours[3*k+j];
491                //    if(useint>0){
492                //    // neighbour is not a boundary condition
493                //        if(stage_centroid_values[useint] > bed_centroid_values[useint]){
494                //            //wet
495                //            neighbours_wet[j]=1;
496                //            useint = useint*3 + neighbour_edges[3*k+j]; // Index of neighbouring edge
497                //            stage_edges[j] = stage_edge_values[useint]; // stage of neighbouring edge
498                //        }else{
499                //            //dry
500                //            neighbours_wet[j]=0;
501                //            stage_edges[j]=0.0; // Arbitrary
502                //        }
503                //    }else{
504                //        // neighbour is a boundary condition, ASSUME IT IS WET
505                //        useint = -useint -1;
506                //        neighbours_wet[j]=1;
507                //        stage_edges[j] = stage_boundary_values[useint];
508                //        }
509                //}
510                ////calculate average neighbouring edge stage, only considering wet triangles
511                //if(max(neighbours_wet[0], max(neighbours_wet[1],neighbours_wet[2]))>0){
512                //    //usedouble == average edge stage
513                //    usedouble = stage_edges[0]*neighbours_wet[0] + stage_edges[1]*neighbours_wet[1] + stage_edges[2]*neighbours_wet[2];
514                //    usedouble /= 1.0*(neighbours_wet[0]+neighbours_wet[1]+neighbours_wet[2]);
515                //}else{
516                //    usedouble = zl;
517                //}
518                //usedouble=min(usedouble,stage_centroid_values[k]);
519                //xmom_explicit_update[k] -= -0.5*g*max(usedouble-zl,0.0)*(usedouble-zl)*normals[ki2]*length;
520                //ymom_explicit_update[k] -= -0.5*g*max(usedouble-zl,0.0)*(usedouble-zl)*normals[ki2+1]*length;
521                ////xmom_explicit_update[k] -= -0.5*g*abs(usedouble-zl)*(usedouble-zl)*normals[ki2]*length;
522                ////ymom_explicit_update[k] -= -0.5*g*abs(usedouble-zl)*(usedouble-zl)*normals[ki2+1]*length;
523
524            //}
525
526            already_computed_flux[ki] = call; // #k Done
527
528
529            // Update neighbour n with same flux but reversed sign
530            if (n >= 0) {
531                stage_explicit_update[n] += edgeflux[0];
532                xmom_explicit_update[n] += edgeflux[1];
533                ymom_explicit_update[n] += edgeflux[2];
534                //Add bed slope term here
535                if((stage_centroid_values[n]>bed_centroid_values[n]+0.0)){
536                    bedslope_work = -0.5*g*max(stage_centroid_values[n]-zr,0.)*(stage_centroid_values[n]-zr)*length;
537                    xmom_explicit_update[n] += normals[ki2]*bedslope_work;
538                    ymom_explicit_update[n] += normals[ki2+1]*bedslope_work;
539                }else{
540                    // Treat nearly dry cells
541                    //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;
542                    bedslope_work = -pressure_flux*length;
543                    xmom_explicit_update[n] += normals[ki2]*bedslope_work;
544                    ymom_explicit_update[n] += normals[ki2+1]*bedslope_work;
545                }
546                    // The triangle is (partially?) dry
547                    // Find the average stage of the wet neighbours edges -- set this as the within triangle stage for the bedslope term
548                    // FIXME: Note that we keep recomputing this -- could improve efficiency
549                    //for(j=0;j<3;j++){
550                    //    // For each edge, check if each neighbouring triangle is wet
551                    //    useint = neighbours[3*n+j];
552                    //    if(useint>0){
553                    //    // neighbour is not a boundary condition
554                    //        if(stage_centroid_values[useint] > bed_centroid_values[useint]){
555                    //            //wet
556                    //            neighbours_wet[j]=1;
557                    //            useint = useint*3 + neighbour_edges[3*n+j]; // Index of neighbouring edge
558                    //            stage_edges[j] = stage_edge_values[useint]; // stage of neighbouring edge
559                    //        }else{
560                    //            //dry
561                    //            neighbours_wet[j]=0;
562                    //            stage_edges[j]=0.0; //Arbitrary
563                    //        }
564                    //    }else{
565                    //        // neighbour is a boundary condition, ASSUME IT IS WET
566                    //        useint = -useint -1;
567                    //        neighbours_wet[j]=1;
568                    //        stage_edges[j] = stage_boundary_values[useint];
569                    //        }
570                    //}
571                    ////calculate average neighbouring edge stage, only considering wet triangles
572                    //if(max(neighbours_wet[0],max(neighbours_wet[1],neighbours_wet[2]))>0){
573                    //    //average edge stage
574                    //    usedouble = stage_edges[0]*neighbours_wet[0] + stage_edges[1]*neighbours_wet[1] + stage_edges[2]*neighbours_wet[2];
575                    //    usedouble /= 1.0*(neighbours_wet[0]+neighbours_wet[1]+neighbours_wet[2]);
576                    //}else{
577                    //    usedouble = zr;
578                    //}
579                    //usedouble=min(usedouble,stage_centroid_values[n]);
580                    //xmom_explicit_update[n] += -0.5*g*max(usedouble-zr,0.0)*(usedouble-zr)*normals[ki2]*length;
581                    //ymom_explicit_update[n] += -0.5*g*max(usedouble-zr,0.0)*(usedouble-zr)*normals[ki2+1]*length;
582                    ////xmom_explicit_update[n] += -0.5*g*abs(usedouble-zr)*(usedouble-zr)*normals[ki2]*length;
583                    ////ymom_explicit_update[n] += -0.5*g*abs(usedouble-zr)*(usedouble-zr)*normals[ki2+1]*length;
584
585               //}
586
587                already_computed_flux[nm] = call; // #n Done
588            }
589
590            // Update timestep based on edge i and possibly neighbour n
591            if (tri_full_flag[k] == 1) {
592                if (max_speed > epsilon) {
593                    // Apply CFL condition for triangles joining this edge (triangle k and triangle n)
594
595                    // CFL for triangle k
596                    timestep = min(timestep, radii[k] / max_speed);
597
598                    if (n >= 0) {
599                        // Apply CFL condition for neigbour n (which is on the ith edge of triangle k)
600                        timestep = min(timestep, radii[n] / max_speed);
601                    }
602
603                    // Ted Rigby's suggested less conservative version
604                    //if(n>=0){
606                    //}else{
607                    //  timestep = min(timestep, radii[k]/max_speed);
608                    //}
609                }
610            }
611
612        } // End edge i (and neighbour n)
613
614
615        // Normalise triangle k by area and store for when all conserved
616        // quantities get updated
617        inv_area = 1.0 / areas[k];
618        stage_explicit_update[k] *= inv_area;
619        xmom_explicit_update[k] *= inv_area;
620        ymom_explicit_update[k] *= inv_area;
621
622
623        // Keep track of maximal speeds
624        max_speed_array[k] = max_speed;
625
626    } // End triangle k
627
628
629    // Try to adjust stage_explicit_update to account for the negative stages,
630    // using an approach assumes a very small discharge is enough to allow
631    // the depth to increase from negative to zero
632    //for(k=0;k<number_of_elements;k++){
633        //depthtemp = stage_centroid_values[k]-bed_centroid_values[k];
634    //    depthtemp = stage_centroid_values[k]-max(stage_edge_values[3*k+1], max(stage_edge_values[3*k+2], stage_edge_values[3*k]));
635    //   if(depthtemp< 0.0){
636    //       depthtemp = -depthtemp; // stage defecit
637    //       usedouble = 2.0; //1.0e+2; //Large constant
638    //       if(stage_explicit_update[k]*timestep*usedouble > depthtemp){
639    //            stage_explicit_update[k] += (depthtemp - depthtemp/usedouble)/timestep;
640    //        }else{
641    //            stage_explicit_update[k] *=usedouble ;
642    //        }
643    //    }
644    //}
645
646    return timestep;
647}
648
649// Protect against the water elevation falling below the triangle bed
650int _protect(int N,
651         double minimum_allowed_height,
652         double maximum_allowed_speed,
653         double epsilon,
654         double* wc,
655         double* zc,
656         double* zv,
657         double* xmomc,
658         double* ymomc) {
659
660  int k;
661  double hc, wmin;
662  double u, v, reduced_speed;
663
664
665  // Protect against initesimal and negative heights
666  if (maximum_allowed_speed < epsilon) {
667    for (k=0; k<N; k++) {
668      hc = wc[k] - zc[k];
669
670      if (hc < minimum_allowed_height) {
671
672        // Set momentum to zero and ensure h is non negative
673        xmomc[k] = 0.0;
674        ymomc[k] = 0.0;
675        if (hc <= 0.0){
676             // Definine the minimum possible stage value as the minimum bed edge value on triangle k.
677             wmin = 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])));
678             wc[k] = max(wc[k], wmin);
679        }
680      }
681    }
682  } else {
683
684    // Protect against initesimal and negative heights
685    for (k=0; k<N; k++) {
686      hc = wc[k] - zc[k];
687
688      if (hc < minimum_allowed_height) {
689
690        //New code: Adjust momentum to guarantee speeds are physical
691        //          ensure h is non negative
692
693        if (hc <= 0.0) {
694            wc[k] = zc[k];
695        xmomc[k] = 0.0;
696        ymomc[k] = 0.0;
697        } else {
698          //Reduce excessive speeds derived from division by small hc
699          //FIXME (Ole): This may be unnecessary with new slope limiters
700          //in effect.
701
702          u = xmomc[k]/hc;
703          if (fabs(u) > maximum_allowed_speed) {
704            reduced_speed = maximum_allowed_speed * u/fabs(u);
705            //printf("Speed (u) has been reduced from %.3f to %.3f\n",
706            //   u, reduced_speed);
707            xmomc[k] = reduced_speed * hc;
708          }
709
710          v = ymomc[k]/hc;
711          if (fabs(v) > maximum_allowed_speed) {
712            reduced_speed = maximum_allowed_speed * v/fabs(v);
713            //printf("Speed (v) has been reduced from %.3f to %.3f\n",
714            //   v, reduced_speed);
715            ymomc[k] = reduced_speed * hc;
716          }
717        }
718      }
719    }
720  }
721  return 0;
722}
723
724//=========================================================================
725// Python Glue
726//=========================================================================
727
728
729//========================================================================
730// Compute fluxes
731//========================================================================
732
733// Modified central flux function
734
735PyObject *compute_fluxes_ext_central(PyObject *self, PyObject *args) {
736  /*Compute all fluxes and the timestep suitable for all volumes
737    in domain.
738
739    Compute total flux for each conserved quantity using "flux_function_central"
740
741    Fluxes across each edge are scaled by edgelengths and summed up
742    Resulting flux is then scaled by area and stored in
743    explicit_update for each of the three conserved quantities
744    stage, xmomentum and ymomentum
745
746    The maximal allowable speed computed by the flux_function for each volume
747    is converted to a timestep that must not be exceeded. The minimum of
748    those is computed as the next overall timestep.
749
750    Python call:
751    domain.timestep = compute_fluxes(timestep,
752                                     domain.epsilon,
753                     domain.H0,
754                                     domain.g,
755                                     domain.neighbours,
756                                     domain.neighbour_edges,
757                                     domain.normals,
758                                     domain.edgelengths,
760                                     domain.areas,
761                                     tri_full_flag,
762                                     Stage.edge_values,
763                                     Xmom.edge_values,
764                                     Ymom.edge_values,
765                                     Bed.edge_values,
766                                     Stage.boundary_values,
767                                     Xmom.boundary_values,
768                                     Ymom.boundary_values,
769                                     Stage.explicit_update,
770                                     Xmom.explicit_update,
771                                     Ymom.explicit_update,
773                                     optimise_dry_cells,
774                                     stage.centroid_values,
775                                     bed.centroid_values)
776
777
778    Post conditions:
779      domain.explicit_update is reset to computed flux values
780      domain.timestep is set to the largest step satisfying all volumes.
781
782
783  */
784
785
786  PyArrayObject *neighbours, *neighbour_edges,
788    *tri_full_flag,
789    *stage_edge_values,
790    *xmom_edge_values,
791    *ymom_edge_values,
792    *bed_edge_values,
793    *stage_boundary_values,
794    *xmom_boundary_values,
795    *ymom_boundary_values,
796    *stage_explicit_update,
797    *xmom_explicit_update,
798    *ymom_explicit_update,
799    *already_computed_flux, //Tracks whether the flux across an edge has already been computed
800    *max_speed_array, //Keeps track of max speeds for each triangle
801    *stage_centroid_values,
802    *bed_centroid_values,
803    *bed_vertex_values;
804
805  double timestep, epsilon, H0, g;
806  int optimise_dry_cells;
807
808  // Convert Python arguments to C
809  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOiOOO",
810            &timestep,
811            &epsilon,
812            &H0,
813            &g,
814            &neighbours,
815            &neighbour_edges,
816            &normals,
818            &tri_full_flag,
819            &stage_edge_values,
820            &xmom_edge_values,
821            &ymom_edge_values,
822            &bed_edge_values,
823            &stage_boundary_values,
824            &xmom_boundary_values,
825            &ymom_boundary_values,
826            &stage_explicit_update,
827            &xmom_explicit_update,
828            &ymom_explicit_update,
830            &max_speed_array,
831            &optimise_dry_cells,
832            &stage_centroid_values,
833            &bed_centroid_values,
834            &bed_vertex_values)) {
835    report_python_error(AT, "could not parse input arguments");
836    return NULL;
837  }
838
839  // check that numpy array objects arrays are C contiguous memory
840  CHECK_C_CONTIG(neighbours);
841  CHECK_C_CONTIG(neighbour_edges);
842  CHECK_C_CONTIG(normals);
843  CHECK_C_CONTIG(edgelengths);
845  CHECK_C_CONTIG(areas);
846  CHECK_C_CONTIG(tri_full_flag);
847  CHECK_C_CONTIG(stage_edge_values);
848  CHECK_C_CONTIG(xmom_edge_values);
849  CHECK_C_CONTIG(ymom_edge_values);
850  CHECK_C_CONTIG(bed_edge_values);
851  CHECK_C_CONTIG(stage_boundary_values);
852  CHECK_C_CONTIG(xmom_boundary_values);
853  CHECK_C_CONTIG(ymom_boundary_values);
854  CHECK_C_CONTIG(stage_explicit_update);
855  CHECK_C_CONTIG(xmom_explicit_update);
856  CHECK_C_CONTIG(ymom_explicit_update);
858  CHECK_C_CONTIG(max_speed_array);
859  CHECK_C_CONTIG(stage_centroid_values);
860  CHECK_C_CONTIG(bed_centroid_values);
861  CHECK_C_CONTIG(bed_vertex_values);
862
863  int number_of_elements = stage_edge_values -> dimensions[0];
864
865  // Call underlying flux computation routine and update
866  // the explicit update arrays
867  timestep = _compute_fluxes_central(number_of_elements,
868                     timestep,
869                     epsilon,
870                     H0,
871                     g,
872                     (long*) neighbours -> data,
873                     (long*) neighbour_edges -> data,
874                     (double*) normals -> data,
875                     (double*) edgelengths -> data,
877                     (double*) areas -> data,
878                     (long*) tri_full_flag -> data,
879                     (double*) stage_edge_values -> data,
880                     (double*) xmom_edge_values -> data,
881                     (double*) ymom_edge_values -> data,
882                     (double*) bed_edge_values -> data,
883                     (double*) stage_boundary_values -> data,
884                     (double*) xmom_boundary_values -> data,
885                     (double*) ymom_boundary_values -> data,
886                     (double*) stage_explicit_update -> data,
887                     (double*) xmom_explicit_update -> data,
888                     (double*) ymom_explicit_update -> data,
890                     (double*) max_speed_array -> data,
891                     optimise_dry_cells,
892                     (double*) stage_centroid_values -> data,
893                     (double*) bed_centroid_values -> data,
894                     (double*) bed_vertex_values -> data);
895
896  // Return updated flux timestep
897  return Py_BuildValue("d", timestep);
898}
899
900
901PyObject *flux_function_central(PyObject *self, PyObject *args) {
902  //
903  // Gateway to innermost flux function.
904  // This is only used by the unit tests as the c implementation is
905  // normally called by compute_fluxes in this module.
906
907
908  PyArrayObject *normal, *ql, *qr,  *edgeflux;
909  double g, epsilon, max_speed, H0, zl, zr;
910  double h0, limiting_threshold, pressure_flux;
911
912  if (!PyArg_ParseTuple(args, "OOOddOddd",
913            &normal, &ql, &qr, &zl, &zr, &edgeflux,
914            &epsilon, &g, &H0)) {
915      report_python_error(AT, "could not parse input arguments");
916      return NULL;
917  }
918
919
920  h0 = H0*H0; // This ensures a good balance when h approaches H0.
921              // But evidence suggests that h0 can be as little as
922              // epsilon!
923
924  limiting_threshold = 10*H0; // Avoid applying limiter below this
925                              // threshold for performance reasons.
926                              // See ANUGA manual under flux limiting
927
928  pressure_flux = 0.0; // Store the water pressure related component of the flux
929  _flux_function_central((double*) ql -> data,
930                         (double*) qr -> data,
931                         zl,
932                         zr,
933                         ((double*) normal -> data)[0],
934                         ((double*) normal -> data)[1],
935                         epsilon, h0, limiting_threshold,
936                         g,
937                         (double*) edgeflux -> data,
938                         &max_speed,
939             &pressure_flux);
940
941  return Py_BuildValue("d", max_speed);
942}
943
944//========================================================================
945// Gravity
946//========================================================================
947
948PyObject *gravity(PyObject *self, PyObject *args) {
949  //
950  //  gravity(g, h, v, x, xmom, ymom)
951  //
952
953
954  PyArrayObject *h, *v, *x, *xmom, *ymom;
955  int k, N, k3, k6;
956  double g, avg_h, zx, zy;
957  double x0, y0, x1, y1, x2, y2, z0, z1, z2;
958  //double epsilon;
959
960  if (!PyArg_ParseTuple(args, "dOOOOO",
961                        &g, &h, &v, &x,
962                        &xmom, &ymom)) {
963    //&epsilon)) {
964    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: gravity could not parse input arguments");
965    return NULL;
966  }
967
968  // check that numpy array objects arrays are C contiguous memory
969  CHECK_C_CONTIG(h);
970  CHECK_C_CONTIG(v);
971  CHECK_C_CONTIG(x);
972  CHECK_C_CONTIG(xmom);
973  CHECK_C_CONTIG(ymom);
974
975  N = h -> dimensions[0];
976  for (k=0; k<N; k++) {
977    k3 = 3*k;  // base index
978
979    // Get bathymetry
980    z0 = ((double*) v -> data)[k3 + 0];
981    z1 = ((double*) v -> data)[k3 + 1];
982    z2 = ((double*) v -> data)[k3 + 2];
983
984    // Optimise for flat bed
985    // Note (Ole): This didn't produce measurable speed up.
986    // Revisit later
987    //if (fabs(z0-z1)<epsilon && fabs(z1-z2)<epsilon) {
988    //  continue;
989    //}
990
991    // Get average depth from centroid values
992    avg_h = ((double *) h -> data)[k];
993
994    // Compute bed slope
995    k6 = 6*k;  // base index
996
997    x0 = ((double*) x -> data)[k6 + 0];
998    y0 = ((double*) x -> data)[k6 + 1];
999    x1 = ((double*) x -> data)[k6 + 2];
1000    y1 = ((double*) x -> data)[k6 + 3];
1001    x2 = ((double*) x -> data)[k6 + 4];
1002    y2 = ((double*) x -> data)[k6 + 5];
1003
1004
1005    _gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2, &zx, &zy);
1006
1007    // Update momentum
1008    ((double*) xmom -> data)[k] += -g*zx*avg_h;
1009    ((double*) ymom -> data)[k] += -g*zy*avg_h;
1010  }
1011
1012  return Py_BuildValue("");
1013}
1014
1015
1016//========================================================================
1017// Protect -- to prevent the water level from falling below the minimum
1018// bed_edge_value
1019//========================================================================
1020
1021PyObject *protect(PyObject *self, PyObject *args) {
1022  //
1023  //    protect(minimum_allowed_height, maximum_allowed_speed, wc, zc, xmomc, ymomc)
1024
1025
1026  PyArrayObject
1027  *wc,            //Stage at centroids
1028  *zc,            //Elevation at centroids
1029  *zv,            //Elevation at vertices
1030  *xmomc,         //Momentums at centroids
1031  *ymomc;
1032
1033
1034  int N;
1035  double minimum_allowed_height, maximum_allowed_speed, epsilon;
1036
1037  // Convert Python arguments to C
1038  if (!PyArg_ParseTuple(args, "dddOOOOO",
1039            &minimum_allowed_height,
1040            &maximum_allowed_speed,
1041            &epsilon,
1042            &wc, &zc, &zv, &xmomc, &ymomc)) {
1043    report_python_error(AT, "could not parse input arguments");
1044    return NULL;
1045  }
1046
1047  N = wc -> dimensions[0];
1048
1049  _protect(N,
1050       minimum_allowed_height,
1051       maximum_allowed_speed,
1052       epsilon,
1053       (double*) wc -> data,
1054       (double*) zc -> data,
1055       (double*) zv -> data,
1056       (double*) xmomc -> data,
1057       (double*) ymomc -> data);
1058
1059  return Py_BuildValue("");
1060}
1061
1062//========================================================================
1063// Method table for python module
1064//========================================================================
1065
1066static struct PyMethodDef MethodTable[] = {
1067  /* The cast of the function is necessary since PyCFunction values
1068   * only take two PyObject* parameters, and rotate() takes
1069   * three.
1070   */
1071  {"compute_fluxes_ext_central", compute_fluxes_ext_central, METH_VARARGS, "Print out"},
1072  {"gravity_c",        gravity,            METH_VARARGS, "Print out"},
1073  {"flux_function_central", flux_function_central, METH_VARARGS, "Print out"},
1074  {"protect",          protect, METH_VARARGS | METH_KEYWORDS, "Print out"},
1075  {NULL, NULL, 0, NULL}
1076};
1077
1078// Module initialisation
1079void initswb2_domain_ext(void){
1080  Py_InitModule("swb2_domain_ext", MethodTable);
1081
1082  import_array(); // Necessary for handling of NumPY structures
1083}
Note: See TracBrowser for help on using the repository browser.