source: anuga_work/development/anuga_1d/comp_flux_ext_steve.c @ 7793

Last change on this file since 7793 was 6454, checked in by steve, 15 years ago

Added Padarn's C code
File size: 11.0 KB
RevLine 
[6454]1#include "Python.h"
2#include "Numeric/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//Innermost flux function (using w=z+h)
13int _flux_function(double *q_left, double *q_right,
14                   double normals, double g, double epsilon, double h0,
15                   double *edgeflux, double *max_speed) {
16               
17                int i;
18                double flux_left[2], flux_right[2];
19                double w_left, h_left, uh_left, z_left, u_left, soundspeed_left;
20                double w_right, h_right, uh_right, z_right, u_right, soundspeed_right;
21                double z, s_max, s_min, denom;
22               
23                //printf("h0 = %f \n",h0);
24                w_left  = q_left[0];
25                uh_left = q_left[1]*normals;
26                h_left  = q_left[2];
27                z_left  = q_left[3];
28                u_left  = q_left[4]*normals;
29       
30                w_right  = q_right[0];
31                uh_right = q_right[1]*normals;
32                h_right  = q_right[2];
33                z_right  = q_right[3];
34                u_right  = q_right[4]*normals;         
35         
36                z = (z_left+z_right)/2.0;                 
37                                   
38                soundspeed_left = sqrt(g*h_left);
39                soundspeed_right = sqrt(g*h_right);
40               
41                s_max = max(u_left+soundspeed_left, u_right+soundspeed_right);
42                if (s_max < 0.0) s_max = 0.0;
43       
44                s_min = min(u_left-soundspeed_left, u_right-soundspeed_right);
45                if (s_min > 0.0) s_min = 0.0;
46               
47               
48                // Flux formulas
49                flux_left[0] = u_left*h_left;
50                flux_left[1] = u_left*u_left*h_left + 0.5*g*h_left*h_left;
51
52                flux_right[0] = u_right*h_right;
53                flux_right[1] = u_right*u_right*h_right + 0.5*g*h_right*h_right;
54
55                // Flux computation
56                denom = s_max-s_min;
57                if (denom < epsilon) {
58                        for (i=0; i<2; i++) edgeflux[i] = 0.0;
59                        *max_speed = 0.0;
60                } else {
61                        edgeflux[0] = s_max*flux_left[0] - s_min*flux_right[0];
62                        edgeflux[0] += s_max*s_min*(q_right[0]-q_left[0]);
63                        edgeflux[0] /= denom;
64                        edgeflux[1] = s_max*flux_left[1] - s_min*flux_right[1];
65                        edgeflux[1] += s_max*s_min*(q_right[1]-q_left[1]);
66                        edgeflux[1] /= denom;
67                        edgeflux[1] *= normals;
68   
69                        // Maximal wavespeed
70                        *max_speed = max(fabs(s_max), fabs(s_min));
71                } 
72                return 0;               
73}
74               
75               
76       
77       
78// Computational function for flux computation
79double _compute_fluxes_ext(
80                           double cfl,
81                           double timestep,
82                           double epsilon,
83                           double g,
84                           double h0,
85                           long* neighbours,
86                           long* neighbour_vertices,
87                           double* normals,
88                           double* areas,
89                           double* stage_edge_values,
90                           double* xmom_edge_values,
91                           double* bed_edge_values,
92                           double* height_edge_values,
93                           double* velocity_edge_values,
94                           double* stage_boundary_values,
95                           double* xmom_boundary_values,
96                           double* bed_boundary_values,
97                           double* height_boundary_values,
98                           double* velocity_boundary_values,
99                           double* stage_explicit_update,
100                           double* xmom_explicit_update,
101                           int number_of_elements,
102                           double* max_speed_array) {
103               
104                double flux[2], ql[5], qr[5], edgeflux[2];
105                double max_speed, normal;
106                int k, i, ki, n, m, nm=0;
107               
108               
109                for (k=0; k<number_of_elements; k++) {
110                        flux[0] = 0.0;
111                        flux[1] = 0.0;
112                       
113                        for (i=0; i<2; i++) {
114                                ki = k*2+i;
115                               
116                                ql[0] = stage_edge_values[ki];
117                                ql[1] = xmom_edge_values[ki];
118                                ql[2] = bed_edge_values[ki];
119                                ql[3] = height_edge_values[ki];
120                                ql[4] = velocity_edge_values[ki];
121                               
122                                n = neighbours[ki];
123                                if (n<0) {
124                                        m = -n-1;
125                                        qr[0] = stage_boundary_values[m];
126                                        qr[1] = xmom_boundary_values[m];
127                                        qr[2] = bed_edge_values[m];
128                                        qr[3] = height_edge_values[m];
129                                        qr[4] = velocity_edge_values[m];
130                                } else {
131                                        m = neighbour_vertices[ki];
132                                        nm = n*2+m;
133                                        qr[0] = stage_edge_values[nm];
134                                        qr[1] = xmom_edge_values[nm];
135                                        qr[2] = bed_edge_values[nm];
136                                        qr[3] = height_edge_values[nm];
137                                        qr[4] = velocity_edge_values[nm];
138                                }
139                               
140                                normal = normals[ki];
141                                _flux_function(ql, qr, normal, g, epsilon, h0, edgeflux, &max_speed);
142                                flux[0] -= edgeflux[0];
143                                flux[1] -= edgeflux[1];
144                               
145                                // Update timestep based on edge i and possibly neighbour n
146                                if (max_speed > epsilon) {
147                                    // Original CFL calculation
148                                   
149                                    timestep = min(timestep, 0.5*cfl*areas[k]/max_speed); 
150                                    if (n>=0) {
151                                        timestep = min(timestep, 0.5*cfl*areas[n]/max_speed); 
152                                    }
153                                }
154                        } // End edge i (and neighbour n)
155                        flux[0] /= areas[k];
156                        stage_explicit_update[k] = flux[0];
157                        flux[1] /= areas[k];
158                        xmom_explicit_update[k] = flux[1];
159                       
160                        //Keep track of maximal speeds
161                        max_speed_array[k]=max_speed;
162                }
163                return timestep;       
164        }
165
166//=========================================================================
167// Python Glue
168//=========================================================================
169PyObject *compute_fluxes_ext(PyObject *self, PyObject *args) {
170 
171    PyArrayObject
172        *neighbours, 
173        *neighbour_vertices,
174        *normals, 
175        *areas,
176        *stage_vertex_values,
177        *xmom_vertex_values,
178        *bed_vertex_values,
179        *height_vertex_values,
180        *velocity_vertex_values,
181        *stage_boundary_values,
182        *xmom_boundary_values,
183        *bed_boundary_values,
184        *height_boundary_values,
185        *velocity_boundary_values,
186        *stage_explicit_update,
187        *xmom_explicit_update,
188        *max_speed_array;
189
190   
191    double timestep, epsilon, g, h0, cfl;
192    int number_of_elements;
193   
194    // Convert Python arguments to C
195    if (!PyArg_ParseTuple(args, "dddddOOOOOOOOOOOOOOOOiO",
196                          &cfl,
197                          &timestep,
198                          &epsilon,
199                          &g,
200                          &h0,
201                          &neighbours,
202                          &neighbour_vertices,
203                          &normals,
204                          &areas,
205                          &stage_vertex_values,
206                          &xmom_vertex_values,
207                          &bed_vertex_values,
208                          &height_vertex_values,
209                          &velocity_vertex_values,
210                          &stage_boundary_values,
211                          &xmom_boundary_values,
212                          &bed_boundary_values,
213                          &height_boundary_values,
214                          &velocity_boundary_values,
215                          &stage_explicit_update,
216                          &xmom_explicit_update,
217                          &number_of_elements,
218                          &max_speed_array)) {
219        PyErr_SetString(PyExc_RuntimeError, "comp_flux_ext.c: compute_fluxes_ext could not parse input");
220        return NULL;
221    }
222
223
224    // Call underlying flux computation routine and update
225    // the explicit update arrays
226    timestep = _compute_fluxes_ext(
227        cfl,
228        timestep,
229        epsilon,
230        g,
231        h0,
232        (long*) neighbours -> data,
233        (long*) neighbour_vertices -> data,
234        (double*) normals -> data,
235        (double*) areas -> data,
236        (double*) stage_vertex_values -> data,
237        (double*) xmom_vertex_values -> data,
238        (double*) bed_vertex_values -> data,
239        (double*) height_vertex_values -> data,
240        (double*) velocity_vertex_values -> data,
241        (double*) stage_boundary_values -> data,
242        (double*) xmom_boundary_values -> data,
243        (double*) bed_boundary_values -> data,
244        (double*) height_boundary_values -> data,
245        (double*) velocity_boundary_values -> data,
246        (double*) stage_explicit_update -> data,
247        (double*) xmom_explicit_update -> data,
248        number_of_elements,
249        (double*) max_speed_array -> data);
250
251
252    // Return updated flux timestep
253    return Py_BuildValue("d", timestep);
254}
255
256
257PyObject *compute_fluxes_ext_short(PyObject *self, PyObject *args) {
258 
259    PyObject
260        *domain,
261        *stage, 
262        *xmom, 
263        *bed,
264        *height,
265        *velocity;
266
267    PyArrayObject
268        *neighbours, 
269        *neighbour_vertices,
270        *normals, 
271        *areas,
272        *stage_vertex_values,
273        *xmom_vertex_values,
274        *bed_vertex_values,
275        *height_vertex_values,
276        *velocity_vertex_values,
277        *stage_boundary_values,
278        *xmom_boundary_values,
279        *bed_boundary_values,
280        *height_boundary_values,
281        *velocity_boundary_values,
282        *stage_explicit_update,
283        *xmom_explicit_update,
284        *max_speed_array;
285   
286  double timestep, epsilon, g, h0, cfl;
287  int number_of_elements;
288
289   
290  // Convert Python arguments to C
291  if (!PyArg_ParseTuple(args, "dOOOOOO",
292                        &timestep,
293                        &domain,
294                        &stage,
295                        &xmom,
296                        &bed,
297                        &height,
298                        &velocity)) {
299      PyErr_SetString(PyExc_RuntimeError, "comp_flux_ext.c: compute_fluxes_ext_short could not parse input");
300      return NULL;
301  }
302
303
304    epsilon           = get_python_double(domain,"epsilon");
305    g                 = get_python_double(domain,"g");
306    h0                = get_python_double(domain,"h0");
307    cfl               = get_python_double(domain,"cfl");
308 
309   
310    neighbours        = get_consecutive_array(domain, "neighbours");
311    neighbour_vertices= get_consecutive_array(domain, "neighbour_vertices"); 
312    normals           = get_consecutive_array(domain, "normals");
313    areas             = get_consecutive_array(domain, "areas");   
314    max_speed_array   = get_consecutive_array(domain, "max_speed_array");
315   
316    stage_vertex_values      = get_consecutive_array(stage,    "vertex_values");   
317    xmom_vertex_values       = get_consecutive_array(xmom,     "vertex_values");   
318    bed_vertex_values        = get_consecutive_array(bed,      "vertex_values");   
319    height_vertex_values     = get_consecutive_array(height,   "vertex_values");   
320    velocity_vertex_values   = get_consecutive_array(velocity, "vertex_values");   
321
322    stage_boundary_values     = get_consecutive_array(stage,     "boundary_values");   
323    xmom_boundary_values      = get_consecutive_array(xmom,      "boundary_values");   
324    bed_boundary_values       = get_consecutive_array(bed,       "boundary_values");   
325    height_boundary_values    = get_consecutive_array(height,    "boundary_values");   
326    velocity_boundary_values  = get_consecutive_array(velocity,  "boundary_values");   
327
328
329    stage_explicit_update = get_consecutive_array(stage, "explicit_update");   
330    xmom_explicit_update  = get_consecutive_array(xmom,  "explicit_update");   
331
332    number_of_elements = stage_vertex_values -> dimensions[0];
333 
334    // Call underlying flux computation routine and update
335    // the explicit update arrays
336    timestep = _compute_fluxes_ext(
337        cfl,
338        timestep,
339        epsilon,
340        g,
341        h0,
342        (long*) neighbours -> data,
343        (long*) neighbour_vertices -> data,
344        (double*) normals -> data,
345        (double*) areas -> data,
346        (double*) stage_vertex_values -> data,
347        (double*) xmom_vertex_values -> data,
348        (double*) bed_vertex_values -> data,
349        (double*) height_vertex_values -> data,
350        (double*) velocity_vertex_values -> data,
351        (double*) stage_boundary_values -> data,
352        (double*) xmom_boundary_values -> data,
353        (double*) bed_boundary_values -> data,
354        (double*) height_boundary_values -> data,
355        (double*) velocity_boundary_values -> data,
356        (double*) stage_explicit_update -> data,
357        (double*) xmom_explicit_update -> data,
358        number_of_elements,
359        (double*) max_speed_array -> data);
360   
361
362  Py_DECREF(neighbours);
363  Py_DECREF(neighbour_vertices);
364  Py_DECREF(normals);
365  Py_DECREF(areas);
366  Py_DECREF(stage_vertex_values);
367  Py_DECREF(xmom_vertex_values);
368  Py_DECREF(bed_vertex_values);
369  Py_DECREF(height_vertex_values);
370  Py_DECREF(velocity_vertex_values);
371  Py_DECREF(stage_boundary_values);
372  Py_DECREF(xmom_boundary_values);
373  Py_DECREF(bed_boundary_values);
374  Py_DECREF(height_boundary_values);
375  Py_DECREF(velocity_boundary_values);
376  Py_DECREF(stage_explicit_update);
377  Py_DECREF(xmom_explicit_update);
378  Py_DECREF(max_speed_array);
379
380
381  // Return updated flux timestep
382  return Py_BuildValue("d", timestep);
383}
384
385//-------------------------------
386// Method table for python module
387//-------------------------------
388
389static struct PyMethodDef MethodTable[] = {
390  {"compute_fluxes_ext", compute_fluxes_ext, METH_VARARGS, "Print out"},
391  {"compute_fluxes_ext_short", compute_fluxes_ext_short, METH_VARARGS, "Print out"},
392  {NULL, NULL}
393};
394
395// Module initialisation
396void initcomp_flux_ext_steve(void){
397  Py_InitModule("comp_flux_ext_steve", MethodTable);
398  import_array();
399}
Note: See TracBrowser for help on using the repository browser.