source: anuga_work/development/anuga_1d/comp_flux_ext.c @ 5563

Last change on this file since 5563 was 5563, checked in by steve, 16 years ago

Working version of comp_flux_ext

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