source: anuga_work/development/flow_1d/sww_flow/sww_vel_comp_flux_ext.c @ 7830

Last change on this file since 7830 was 7830, checked in by steve, 12 years ago

Moving channel code to numpy

File size: 17.7 KB
Line 
1#include "Python.h"
2#include "numpy/arrayobject.h"
3#include "math.h"
4#include <stdio.h>
5const double pi = 3.14159265358979;
6
7
8// Shared code snippets
9#include "util_ext.h"
10
11
12/* double max(double a, double b) { */
13/*      double z; */
14/*      z=(a>b)?a:b; */
15/*      return z;} */
16       
17/* double min(double a, double b) { */
18/*      double z; */
19/*      z=(a<b)?a:b; */
20/*      return z;} */
21
22
23// Function to obtain speed from momentum and depth.
24// This is used by flux functions
25// Input parameters uh and h may be modified by this function.
26double _compute_speed(double *uh, 
27                      double *h, 
28                      double epsilon, 
29                      double h0) {
30 
31  double u;
32 
33  if (*h < epsilon) {
34    *h = 0.0;  //Could have been negative
35     u = 0.0;
36  } else {
37    u = *uh/(*h + h0/ *h);   
38  }
39 
40
41  // Adjust momentum to be consistent with speed
42  *uh = u * *h;
43 
44  return u;
45}
46
47//-------------------------------------------------------------
48// New vel based code
49//-------------------------------------------------------------
50//Innermost flux function (using w=z+h)
51int _flux_function_vel(double *q_left, double *q_right,
52                       double normals, double g, double epsilon, double h0,
53                       double *edgeflux, double *max_speed) {
54   
55    int i;
56    double flux_left[2], flux_right[2];
57    double w_left, h_left, uh_left, z_left, u_left, soundspeed_left;
58    double w_right, h_right, uh_right, z_right, u_right, soundspeed_right;
59    double z, s_max, s_min, denom;
60   
61
62    w_left  = q_left[0];
63    uh_left = q_left[1]*normals;
64    z_left  = q_left[2];
65    h_left  = q_left[3];
66    u_left  = q_left[4]*normals;
67
68/*     printf("w_left  = %f \n",w_left); */
69/*     printf("uh_left = %f \n",uh_left); */
70/*     printf("z_left  = %f \n",z_left); */
71/*     printf("h_left  = %f \n",h_left); */
72/*     printf("u_left  = %f \n",u_left); */
73   
74    w_right  = q_right[0];
75    uh_right = q_right[1]*normals;
76    z_right  = q_right[2];
77    h_right  = q_right[3];
78    u_right  = q_right[4]*normals;             
79   
80    z = (z_left+z_right)/2.0;             
81   
82    soundspeed_left = sqrt(g*h_left);
83    soundspeed_right = sqrt(g*h_right);
84   
85    s_max = max(u_left+soundspeed_left, u_right+soundspeed_right);
86    if (s_max < 0.0) s_max = 0.0;
87   
88    s_min = min(u_left-soundspeed_left, u_right-soundspeed_right);
89    if (s_min > 0.0) s_min = 0.0;
90   
91   
92    // Flux formulas
93    flux_left[0] = u_left*h_left;
94    flux_left[1] = u_left*u_left*h_left + 0.5*g*h_left*h_left;
95   
96    flux_right[0] = u_right*h_right;
97    flux_right[1] = u_right*u_right*h_right + 0.5*g*h_right*h_right;
98   
99    // Flux computation
100    denom = s_max-s_min;
101    if (denom < epsilon) {
102        for (i=0; i<2; i++) edgeflux[i] = 0.0;
103        *max_speed = 0.0;
104    } else {
105        edgeflux[0] = s_max*flux_left[0] - s_min*flux_right[0];
106        edgeflux[0] += s_max*s_min*(w_right-w_left);
107        edgeflux[0] /= denom;
108        edgeflux[1] = s_max*flux_left[1] - s_min*flux_right[1];
109        edgeflux[1] += s_max*s_min*(uh_right-uh_left);
110        edgeflux[1] /= denom;
111        edgeflux[1] *= normals;
112       
113        // Maximal wavespeed
114        *max_speed = max(fabs(s_max), fabs(s_min));
115    } 
116    return 0;           
117}
118
119// Computational function for flux computation
120double _compute_fluxes_vel_ext(double cfl,
121                               double timestep,
122                               double epsilon,
123                               double g,
124                               double h0,
125                               long* neighbours,
126                               long* neighbour_vertices,
127                               double* normals,
128                               double* areas,
129                               double* stage_edge_values,
130                               double* xmom_edge_values,
131                               double* bed_edge_values,
132                               double* height_edge_values,
133                               double* velocity_edge_values,
134                               double* stage_boundary_values,
135                               double* xmom_boundary_values,
136                               double* bed_boundary_values,
137                               double* height_boundary_values,
138                               double* velocity_boundary_values,
139                               double* stage_explicit_update,
140                               double* xmom_explicit_update,
141                               int number_of_elements,
142                               double* max_speed_array) {
143               
144    double flux[2], ql[5], qr[5], edgeflux[2];
145    double max_speed, normal;
146    int k, i, ki, n, m, nm=0;
147
148    //printf("Inside _comp\n");   
149   
150    for (k=0; k<number_of_elements; k++) {
151        flux[0] = 0.0;
152        flux[1] = 0.0;
153        //printf("k = %d\n",k);
154
155        for (i=0; i<2; i++) {
156            ki = k*2+i;
157           
158            ql[0] = stage_edge_values[ki];
159            ql[1] = xmom_edge_values[ki];
160            ql[2] = bed_edge_values[ki];
161            ql[3] = height_edge_values[ki];
162            ql[4] = velocity_edge_values[ki];
163           
164            n = neighbours[ki];
165            if (n<0) {
166                m = -n-1;
167                qr[0] = stage_boundary_values[m];
168                qr[1] = xmom_boundary_values[m];
169                qr[2] = bed_boundary_values[m];
170                qr[3] = height_boundary_values[m];
171                qr[4] = velocity_boundary_values[m];
172            } else {
173                m = neighbour_vertices[ki];
174                nm = n*2+m;
175                qr[0] = stage_edge_values[nm];
176                qr[1] = xmom_edge_values[nm];
177                qr[2] = bed_edge_values[nm];
178                qr[3] = height_edge_values[nm];
179                qr[4] = velocity_edge_values[nm];
180            }
181           
182            normal = normals[ki];
183            _flux_function_vel(ql, qr, normal, g, epsilon, h0, edgeflux, &max_speed);
184            flux[0] -= edgeflux[0];
185            flux[1] -= edgeflux[1];
186           
187            // Update timestep based on edge i and possibly neighbour n
188            if (max_speed > epsilon) {
189                // Original CFL calculation
190               
191                timestep = min(timestep, 0.5*cfl*areas[k]/max_speed); 
192                if (n>=0) {
193                    timestep = min(timestep, 0.5*cfl*areas[n]/max_speed); 
194                }
195            }
196        } // End edge i (and neighbour n)
197        flux[0] /= areas[k];
198        stage_explicit_update[k] = flux[0];
199        flux[1] /= areas[k];
200        xmom_explicit_update[k] = flux[1];
201       
202        //Keep track of maximal speeds
203        max_speed_array[k]=max_speed;
204    }
205    return timestep;   
206}
207
208//-------------------------------------------------------------
209// Old code
210//------------------------------------------------------------
211//Innermost flux function (using w=z+h)
212int _flux_function(double *q_left, double *q_right,
213        double z_left, double z_right,
214                   double normals, double g, double epsilon, double h0,
215                double *edgeflux, double *max_speed) {
216               
217                int i;
218                double ql[2], qr[2], flux_left[2], flux_right[2];
219                double z, w_left, h_left, uh_left, soundspeed_left, u_left;
220                double w_right, h_right, uh_right, soundspeed_right, u_right;
221                double s_max, s_min, denom;
222               
223                //printf("h0 = %f \n",h0);
224                ql[0] = q_left[0];
225                ql[1] = q_left[1];
226                ql[1] = ql[1]*normals;
227       
228                qr[0] = q_right[0];
229                qr[1] = q_right[1];
230                qr[1] = qr[1]*normals;
231         
232                z = (z_left+z_right)/2.0;                 
233                                   
234                //w_left = ql[0];
235                //h_left = w_left-z;
236                //uh_left = ql[1];
237               
238
239
240                // Compute speeds in x-direction
241                w_left = ql[0];         
242                h_left = w_left-z;
243                uh_left = ql[1];
244
245                u_left = _compute_speed(&uh_left, &h_left, epsilon, h0);
246
247                w_right = qr[0];
248                h_right = w_right-z;
249                uh_right = qr[1];
250
251                u_right = _compute_speed(&uh_right, &h_right, epsilon, h0);
252 
253                soundspeed_left = sqrt(g*h_left);
254                soundspeed_right = sqrt(g*h_right);
255               
256                s_max = max(u_left+soundspeed_left, u_right+soundspeed_right);
257                if (s_max < 0.0) s_max = 0.0;
258       
259                s_min = min(u_left-soundspeed_left, u_right-soundspeed_right);
260                if (s_min > 0.0) s_min = 0.0;
261               
262               
263                // Flux formulas
264                flux_left[0] = u_left*h_left;
265                flux_left[1] = u_left*uh_left + 0.5*g*h_left*h_left;
266
267                flux_right[0] = u_right*h_right;
268                flux_right[1] = u_right*uh_right + 0.5*g*h_right*h_right;
269
270                // Flux computation
271                denom = s_max-s_min;
272                if (denom < epsilon) {
273                        for (i=0; i<2; i++) edgeflux[i] = 0.0;
274                        *max_speed = 0.0;
275                } else {
276                        edgeflux[0] = s_max*flux_left[0] - s_min*flux_right[0];
277                        edgeflux[0] += s_max*s_min*(qr[0]-ql[0]);
278                        edgeflux[0] /= denom;
279                        edgeflux[1] = s_max*flux_left[1] - s_min*flux_right[1];
280                        edgeflux[1] += s_max*s_min*(qr[1]-ql[1]);
281                        edgeflux[1] /= denom;
282                        edgeflux[1] *= normals;
283   
284                // Maximal wavespeed
285        *max_speed = max(fabs(s_max), fabs(s_min));
286                } 
287    return 0;           
288        }
289               
290               
291       
292       
293// Computational function for flux computation
294double _compute_fluxes_ext(
295                           double cfl,
296                           double timestep,
297                           double epsilon,
298                           double g,
299                           double h0,
300                           long* neighbours,
301                           long* neighbour_vertices,
302                           double* normals,
303                           double* areas,
304                           double* stage_edge_values,
305                           double* xmom_edge_values,
306                           double* bed_edge_values,
307                           double* stage_boundary_values,
308                           double* xmom_boundary_values,
309                           double* stage_explicit_update,
310                           double* xmom_explicit_update,
311                           int number_of_elements,
312                           double* max_speed_array) {
313               
314                double flux[2], ql[2], qr[2], edgeflux[2];
315                double zl, zr, max_speed, normal;
316                int k, i, ki, n, m, nm=0;
317               
318               
319                for (k=0; k<number_of_elements; k++) {
320                        flux[0] = 0.0;
321                        flux[1] = 0.0;
322                       
323                        for (i=0; i<2; i++) {
324                                ki = k*2+i;
325                               
326                                ql[0] = stage_edge_values[ki];
327                                ql[1] = xmom_edge_values[ki];
328                                zl = bed_edge_values[ki];
329                               
330                                n = neighbours[ki];
331                                if (n<0) {
332                                        m = -n-1;
333                                        qr[0] = stage_boundary_values[m];
334                                        qr[1] = xmom_boundary_values[m];
335                                        zr = zl;
336                                } else {
337                                        m = neighbour_vertices[ki];
338                                        nm = n*2+m;
339                                        qr[0] = stage_edge_values[nm];
340                                        qr[1] = xmom_edge_values[nm];
341                                        zr = bed_edge_values[nm];                               
342                                }
343                               
344                                normal = normals[ki];
345                                _flux_function(ql, qr, zl, zr, normal, g, epsilon, h0, edgeflux, &max_speed);
346                                flux[0] -= edgeflux[0];
347                                flux[1] -= edgeflux[1];
348                               
349                                // Update timestep based on edge i and possibly neighbour n
350                                if (max_speed > epsilon) {
351                                    // Original CFL calculation
352                                   
353                                    timestep = min(timestep, 0.5*cfl*areas[k]/max_speed); 
354                                    if (n>=0) {
355                                        timestep = min(timestep, 0.5*cfl*areas[n]/max_speed); 
356                                    }
357                                }
358            } // End edge i (and neighbour n)
359                        flux[0] /= areas[k];
360                        stage_explicit_update[k] = flux[0];
361                        flux[1] /= areas[k];
362                        xmom_explicit_update[k] = flux[1];
363                       
364                        //Keep track of maximal speeds
365                        max_speed_array[k]=max_speed;
366                }
367                return timestep;       
368        }
369
370//=========================================================================
371// Python Glue
372//=========================================================================
373PyObject *compute_fluxes_ext(PyObject *self, PyObject *args) {
374 
375   PyObject
376        *domain,
377        *stage, 
378        *xmom, 
379        *bed;
380
381    PyArrayObject
382        *neighbours, 
383        *neighbour_vertices,
384        *normals, 
385        *areas,
386        *stage_vertex_values,
387        *xmom_vertex_values,
388        *bed_vertex_values,
389        *stage_boundary_values,
390        *xmom_boundary_values,
391        *stage_explicit_update,
392        *xmom_explicit_update,
393        *max_speed_array;
394   
395  double timestep, epsilon, g, h0, cfl;
396  int number_of_elements;
397
398   
399  // Convert Python arguments to C
400  if (!PyArg_ParseTuple(args, "dOOOO",
401                        &timestep,
402                        &domain,
403                        &stage,
404                        &xmom,
405                        &bed)) {
406      PyErr_SetString(PyExc_RuntimeError, "comp_flux_vel_ext.c: compute_fluxes_ext could not parse input");
407      return NULL;
408  }
409
410
411    epsilon           = get_python_double(domain,"epsilon");
412    g                 = get_python_double(domain,"g");
413    h0                = get_python_double(domain,"h0");
414    cfl               = get_python_double(domain,"CFL");
415 
416   
417    neighbours        = get_consecutive_array(domain, "neighbours");
418    neighbour_vertices= get_consecutive_array(domain, "neighbour_vertices"); 
419    normals           = get_consecutive_array(domain, "normals");
420    areas             = get_consecutive_array(domain, "areas");   
421    max_speed_array   = get_consecutive_array(domain, "max_speed_array");
422   
423    stage_vertex_values = get_consecutive_array(stage, "vertex_values");   
424    xmom_vertex_values  = get_consecutive_array(xmom, "vertex_values");   
425    bed_vertex_values   = get_consecutive_array(bed, "vertex_values");   
426
427    stage_boundary_values = get_consecutive_array(stage, "boundary_values");   
428    xmom_boundary_values  = get_consecutive_array(xmom, "boundary_values");   
429
430
431    stage_explicit_update = get_consecutive_array(stage, "explicit_update");   
432    xmom_explicit_update  = get_consecutive_array(xmom, "explicit_update");   
433
434
435
436    number_of_elements = stage_vertex_values -> dimensions[0];
437
438
439 
440    // Call underlying flux computation routine and update
441    // the explicit update arrays
442    timestep = _compute_fluxes_ext(
443        cfl,
444        timestep,
445        epsilon,
446        g,
447        h0,
448        (long*) neighbours -> data,
449        (long*) neighbour_vertices -> data,
450        (double*) normals -> data,
451        (double*) areas -> data,
452        (double*) stage_vertex_values -> data,
453        (double*) xmom_vertex_values -> data,
454        (double*) bed_vertex_values -> data,
455        (double*) stage_boundary_values -> data,
456        (double*) xmom_boundary_values -> data,
457        (double*) stage_explicit_update -> data,
458        (double*) xmom_explicit_update -> data,
459        number_of_elements,
460        (double*) max_speed_array -> data);
461
462
463  Py_DECREF(neighbours);
464  Py_DECREF(neighbour_vertices);
465  Py_DECREF(normals);
466  Py_DECREF(areas);
467  Py_DECREF(stage_vertex_values);
468  Py_DECREF(xmom_vertex_values);
469  Py_DECREF(bed_vertex_values);
470  Py_DECREF(stage_boundary_values);
471  Py_DECREF(xmom_boundary_values);
472  Py_DECREF(stage_explicit_update);
473  Py_DECREF(xmom_explicit_update);
474  Py_DECREF(max_speed_array);
475
476
477
478
479  // Return updated flux timestep
480  return Py_BuildValue("d", timestep);
481}
482
483
484//------------------------------------------------
485// New velocity based compute fluxes
486//------------------------------------------------
487PyObject *compute_fluxes_vel_ext(PyObject *self, PyObject *args) {
488 
489    PyObject
490        *domain,
491        *stage, 
492        *xmom, 
493        *bed,
494        *height,
495        *velocity;
496
497    PyArrayObject
498        *neighbours, 
499        *neighbour_vertices,
500        *normals, 
501        *areas,
502        *stage_vertex_values,
503        *xmom_vertex_values,
504        *bed_vertex_values,
505        *height_vertex_values,
506        *velocity_vertex_values,
507        *stage_boundary_values,
508        *xmom_boundary_values,
509        *bed_boundary_values,
510        *height_boundary_values,
511        *velocity_boundary_values,
512        *stage_explicit_update,
513        *xmom_explicit_update,
514        *max_speed_array;
515   
516  double timestep, epsilon, g, h0, cfl;
517  int number_of_elements;
518
519   
520  // Convert Python arguments to C
521  if (!PyArg_ParseTuple(args, "dOOOOOO",
522                        &timestep,
523                        &domain,
524                        &stage,
525                        &xmom,
526                        &bed,
527                        &height,
528                        &velocity)) {
529      PyErr_SetString(PyExc_RuntimeError, "comp_flux_vel_ext.c: compute_fluxes_vel_ext could not parse input");
530      printf("comp_flux_vel_ext.c: compute_fluxes_vel_ext could not parse input");
531      return NULL;
532  }
533
534
535    epsilon           = get_python_double(domain,"epsilon");
536    g                 = get_python_double(domain,"g");
537    h0                = get_python_double(domain,"h0");
538    cfl               = get_python_double(domain,"CFL");
539 
540   
541    neighbours        = get_consecutive_array(domain, "neighbours");
542    neighbour_vertices= get_consecutive_array(domain, "neighbour_vertices"); 
543    normals           = get_consecutive_array(domain, "normals");
544    areas             = get_consecutive_array(domain, "areas");   
545    max_speed_array   = get_consecutive_array(domain, "max_speed_array");
546   
547    stage_vertex_values      = get_consecutive_array(stage,    "vertex_values");   
548    xmom_vertex_values       = get_consecutive_array(xmom,     "vertex_values");   
549    bed_vertex_values        = get_consecutive_array(bed,      "vertex_values");   
550    height_vertex_values     = get_consecutive_array(height,   "vertex_values");   
551    velocity_vertex_values   = get_consecutive_array(velocity, "vertex_values");   
552
553    stage_boundary_values     = get_consecutive_array(stage,     "boundary_values");   
554    xmom_boundary_values      = get_consecutive_array(xmom,      "boundary_values");   
555    bed_boundary_values       = get_consecutive_array(bed,       "boundary_values");   
556    height_boundary_values    = get_consecutive_array(height,    "boundary_values");   
557    velocity_boundary_values  = get_consecutive_array(velocity,  "boundary_values");   
558
559
560    stage_explicit_update = get_consecutive_array(stage, "explicit_update");   
561    xmom_explicit_update  = get_consecutive_array(xmom,  "explicit_update");   
562
563    number_of_elements = stage_vertex_values -> dimensions[0];
564 
565
566    //printf("comp_flux_vel 1 N = %d %15.5e %15.5e %15.5e %15.5e %15.5e \n",number_of_elements,cfl,timestep,epsilon,g,h0);
567
568
569
570    // Call underlying flux computation routine and update
571    // the explicit update arrays
572    timestep = _compute_fluxes_vel_ext(cfl,
573                                       timestep,
574                                       epsilon,
575                                       g,
576                                       h0,
577                                       (long*) neighbours -> data,
578                                       (long*) neighbour_vertices -> data,
579                                       (double*) normals -> data,
580                                       (double*) areas -> data,
581                                       (double*) stage_vertex_values -> data,
582                                       (double*) xmom_vertex_values -> data,
583                                       (double*) bed_vertex_values -> data,
584                                       (double*) height_vertex_values -> data,
585                                       (double*) velocity_vertex_values -> data,
586                                       (double*) stage_boundary_values -> data,
587                                       (double*) xmom_boundary_values -> data,
588                                       (double*) bed_boundary_values -> data,
589                                       (double*) height_boundary_values -> data,
590                                       (double*) velocity_boundary_values -> data,
591                                       (double*) stage_explicit_update -> data,
592                                       (double*) xmom_explicit_update -> data,
593                                       number_of_elements,
594                                       (double*) max_speed_array -> data);
595
596    //printf("comp_flux_vel 2 \n");   
597   
598    Py_DECREF(neighbours);
599    Py_DECREF(neighbour_vertices);
600    Py_DECREF(normals);
601    Py_DECREF(areas);
602    Py_DECREF(stage_vertex_values);
603    Py_DECREF(xmom_vertex_values);
604    Py_DECREF(bed_vertex_values);
605    Py_DECREF(height_vertex_values);
606    Py_DECREF(velocity_vertex_values);
607    Py_DECREF(stage_boundary_values);
608    Py_DECREF(xmom_boundary_values);
609    Py_DECREF(bed_boundary_values);
610    Py_DECREF(height_boundary_values);
611    Py_DECREF(velocity_boundary_values);
612    Py_DECREF(stage_explicit_update);
613    Py_DECREF(xmom_explicit_update);
614    Py_DECREF(max_speed_array);
615   
616    //printf("comp_flux_vel 3 \n");
617
618    // Return updated flux timestep
619    return Py_BuildValue("d", timestep);
620}
621
622
623
624//-------------------------------
625// Method table for python module
626//-------------------------------
627
628static struct PyMethodDef MethodTable[] = {
629  {"compute_fluxes_ext", compute_fluxes_ext, METH_VARARGS, "Print out"},
630  {"compute_fluxes_vel_ext", compute_fluxes_vel_ext, METH_VARARGS, "Print out"},
631  {NULL, NULL}
632};
633
634// Module initialisation
635void initsww_vel_comp_flux_ext(void){
636  Py_InitModule("sww_vel_comp_flux_ext", MethodTable);
637  import_array();
638}
Note: See TracBrowser for help on using the repository browser.