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

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

Adding C extension for flux calculation
File size: 6.7 KB
RevLine 
[5562]1#include "Python.h"
2#include "Numeric/arrayobject.h"
3#include "math.h"
4#include <stdio.h>
5const double pi = 3.14159265358979;
6
7double max(double a, double b) {
8        double z;
9        z=(a>b)?a:b;
10        return z;}
11       
12double min(double a, double b) {
13        double z;
14        z=(a<b)?a:b;
15        return z;}
16       
17
18//Innermost flux function (using w=z+h)
19int _flux_function(double *q_left, double *q_right,
20        double z_left, double z_right,
21                double normals, double g, double epsilon,
22                double *edgeflux, double *max_speed) {
23               
24                int i;
25                double ql[2], qr[2], flux_left[2], flux_right[2];
26                double z, w_left, h_left, uh_left, soundspeed_left, u_left;
27                double w_right, h_right, uh_right, soundspeed_right, u_right;
28                double s_max, s_min, denom;
29               
30               
31                ql[0] = q_left[0];
32                ql[1] = q_left[1];
33                ql[1] = ql[1]*normals;
34       
35                qr[0] = q_right[0];
36                qr[1] = q_right[1];
37                qr[1] = qr[1]*normals;
38         
39                z = (z_left+z_right)/2.0;                 
40                                   
41                w_left = ql[0];
42                h_left = w_left-z;
43                uh_left = ql[1];
44               
45
46
47                // Compute speeds in x-direction
48                w_left = ql[0];         
49                h_left = w_left-z;
50                uh_left = ql[1];
51                if (h_left < epsilon) {
52                        u_left = 0.0; h_left = 0.0;
53                }
54                else {
55                        u_left = uh_left/h_left;
56                }
57       
58                w_right = qr[0];
59                h_right = w_right-z;
60                uh_right = qr[1];
61                if (h_right < epsilon) {
62                        u_right = 0.0; h_right = 0.0; 
63                }
64                else {
65                        u_right = uh_right/h_right;
66                }
67 
68                soundspeed_left = sqrt(g*h_left);
69                soundspeed_right = sqrt(g*h_right);
70               
71                s_max = max(u_left+soundspeed_left, u_right+soundspeed_right);
72                if (s_max < 0.0) s_max = 0.0;
73       
74                s_min = min(u_left-soundspeed_left, u_right-soundspeed_right);
75                if (s_min > 0.0) s_min = 0.0;
76               
77               
78                // Flux formulas
79                flux_left[0] = u_left*h_left;
80                flux_left[1] = u_left*uh_left + 0.5*g*h_left*h_left;
81
82                flux_right[0] = u_right*h_right;
83                flux_right[1] = u_right*uh_right + 0.5*g*h_right*h_right;
84
85                // Flux computation
86                denom = s_max-s_min;
87                if (denom < epsilon) {
88                        for (i=0; i<2; i++) edgeflux[i] = 0.0;
89                        *max_speed = 0.0;
90                } else {
91                        edgeflux[0] = s_max*flux_left[0] - s_min*flux_right[0];
92                        edgeflux[0] += s_max*s_min*(qr[0]-ql[0]);
93                        edgeflux[0] /= denom;
94                        edgeflux[1] = s_max*flux_left[1] - s_min*flux_right[1];
95                        edgeflux[1] += s_max*s_min*(qr[1]-ql[1]);
96                        edgeflux[1] /= denom;
97                        edgeflux[1] *= normals;
98   
99                // Maximal wavespeed
100        *max_speed = max(fabs(s_max), fabs(s_min));
101                } 
102    return 0;           
103        }
104               
105               
106       
107       
108// Computational function for flux computation
109double _compute_fluxes_ext(double timestep,
110                double epsilon,
111                double g,
112                long* neighbours,
113                long* neighbour_vertices,
114                double* normals,
115                double* areas,
116                double* stage_edge_values,
117                double* xmom_edge_values,
118                double* bed_edge_values,
119                double* stage_boundary_values,
120                double* xmom_boundary_values,
121                double* stage_explicit_update,
122                double* xmom_explicit_update,
123                int number_of_elements,
124                double* max_speed_array) {
125               
126                double flux[2], ql[2], qr[2], edgeflux[2];
127                double zl, zr, max_speed, normal;
128                int k, i, ki, n, m, nm=0;
129               
130                for (k=0; k<number_of_elements; k++) {
131                        flux[0] = 0.0;
132                        flux[1] = 0.0;
133                       
134                        for (i=0; i<2; i++) {
135                                ki = k*2+i;
136                               
137                                ql[0] = stage_edge_values[ki];
138                                ql[1] = xmom_edge_values[ki];
139                                zl = bed_edge_values[ki];
140                               
141                                n = neighbours[ki];
142                                if (n<0) {
143                                        m = -n-1;
144                                        qr[0] = stage_boundary_values[m];
145                                        qr[1] = xmom_boundary_values[m];
146                                        zr = zl;
147                                } else {
148                                        m = neighbour_vertices[ki];
149                                        nm = n*2+m;
150                                        qr[0] = stage_edge_values[nm];
151                                        qr[1] = xmom_edge_values[nm];
152                                        zr = bed_edge_values[nm];                               
153                                }
154                               
155                                normal = normals[ki];
156                                _flux_function(ql, qr, zl, zr, normal, g, epsilon, flux, &max_speed);
157                                flux[0] -= edgeflux[0];
158                                flux[1] -= edgeflux[1];
159                               
160                                // Update timestep based on edge i and possibly neighbour n
161                                if (max_speed > epsilon) {
162                        // Original CFL calculation
163                                       
164                                        timestep = min(timestep, 0.5*areas[k]/max_speed); //Here, CFL=1.0 is assumed. ?????????????????????????????????????????????
165                                        if (n>=0) {
166                                                timestep = min(timestep, 0.5*areas[n]/max_speed); //Here, CFL=1.0 is assumed. ?????????????????????????????????????????????
167                                        }
168                                }
169            } // End edge i (and neighbour n)
170                        flux[0] /= areas[k];
171                        stage_explicit_update[k] = flux[0];
172                        flux[1] /= areas[k];
173                        xmom_explicit_update[k] = flux[1];
174                       
175                        //Keep track of maximal speeds
176                        max_speed_array[k]=max_speed;
177                }
178                return timestep;        }
179       
180       
181       
182       
183       
184
185
186
187//=========================================================================
188// Python Glue
189//=========================================================================
190PyObject *compute_fluxes_ext(PyObject *self, PyObject *args) {
191  PyArrayObject *neighbours, 
192        *neighbour_vertices,
193    *normals, 
194        *areas,
195    *stage_edge_values,
196    *xmom_edge_values,
197    *bed_edge_values,
198    *stage_boundary_values,
199    *xmom_boundary_values,
200        *stage_explicit_update,
201        *xmom_explicit_update,
202        *max_speed_array;
203   
204  double timestep, epsilon, g;
205  int number_of_elements;
206   
207  // Convert Python arguments to C
208  if (!PyArg_ParseTuple(args, "dddOOOOOOOOOOOiO",
209                        &timestep,
210                        &epsilon,
211                        &g,
212                        &neighbours,
213                        &neighbour_vertices,
214                        &normals,
215                        &areas,
216                        &stage_edge_values,
217                        &xmom_edge_values,
218                        &bed_edge_values,
219                        &stage_boundary_values,
220                        &xmom_boundary_values,
221                        &stage_explicit_update,
222                        &xmom_explicit_update,
223                        &number_of_elements,
224                        &max_speed_array)) {
225    PyErr_SetString(PyExc_RuntimeError, "comp_flux_ext.c: compute_fluxes_ext could not parse input");
226    return NULL;
227  }
228
229 
230  // Call underlying flux computation routine and update
231  // the explicit update arrays
232  timestep = _compute_fluxes_ext(timestep,
233                                     epsilon,
234                                     g,
235                                     (long*) neighbours -> data,
236                                     (long*) neighbour_vertices -> data,
237                                     (double*) normals -> data,
238                                     (double*) areas -> data,
239                                     (double*) stage_edge_values -> data,
240                                     (double*) xmom_edge_values -> data,
241                                     (double*) bed_edge_values -> data,
242                                     (double*) stage_boundary_values -> data,
243                                     (double*) xmom_boundary_values -> data,
244                                         (double*) stage_explicit_update -> data,
245                                         (double*) xmom_explicit_update -> data,
246                                         number_of_elements,
247                                         (double*) max_speed_array -> data);
248  // Return updated flux timestep
249  return Py_BuildValue("d", timestep);
250}
251
252 
253
254
255//-------------------------------
256// Method table for python module
257//-------------------------------
258
259static struct PyMethodDef MethodTable[] = {
260  {"compute_fluxes_ext", compute_fluxes_ext, METH_VARARGS, "Print out"},
261  {NULL, NULL}
262};
263
264// Module initialisation
265void initcomp_flux_ext(void){
266  Py_InitModule("comp_flux_ext", MethodTable);
267  import_array();
268}
Note: See TracBrowser for help on using the repository browser.