source: anuga_work/development/2010-projects/anuga_1d/sww/swb_comp_flux_ext.c @ 7840

Last change on this file since 7840 was 7839, checked in by steve, 14 years ago

Changing name of 1d projects so that it will be easy to moveto the trunk
File size: 8.8 KB
Line 
1#include "Python.h"
2#include "numpy/arrayobject.h"
3#include "math.h"
4#include <stdio.h>
5#include "util_ext.h"
6const double pi = 3.14159265358979;
7
8//Innermost flux function (using w=z+h)
9int _flux_function(double *q_left, 
10                                                                double *q_right,
11                                                                double normals, 
12                                                                double g, 
13                                                                double epsilon,
14                                                                double *edgeflux, 
15                                                                double *max_speed) {           
16                int i;
17                double z_left, z_right;
18                double flux_left[2], flux_right[2];
19                double z_star, w_left, h_left, h_left_star, uh_left, soundspeed_left, u_left;
20                double w_right, h_right, h_right_star, uh_right, soundspeed_right, u_right;
21                double s_max, s_min, denom;
22               
23                w_left  = q_left[0];
24                uh_left = q_left[1];
25                uh_left = uh_left*normals;
26                z_left  = q_left[2];
27                h_left  = q_left[3];
28                u_left  = q_left[4];
29       
30                w_right  = q_right[0];
31                uh_right = q_right[1];
32                uh_right = uh_right*normals;
33                z_right  = q_right[2];
34                h_right  = q_right[3];
35                u_right  = q_right[4];
36         
37                if (z_left-z_right != 0.0 ) printf("z_l - z_r = %f \n",z_left-z_right);
38       
39                z_star = max(z_left, z_right);                                                                                                          //equation(7) 
40                                   
41                // Compute speeds in x-direction.
42                h_left_star = max(0.0, w_left-z_star);                                                                                                          //equation(8)
43               
44                if (h_left_star < epsilon) {
45                        u_left = 0.0; h_left_star = 0.0;
46                }
47                else {
48                        u_left = uh_left/h_left_star;
49                }
50                                                                                                                                                                                                   
51                h_right_star = max(0.0, w_right-z_star);                                                                                                        //equation(9)
52
53                if (h_right_star < epsilon) {
54                        u_right = 0.0; h_right_star = 0.0; 
55                }
56                else {
57                        u_right = uh_right/h_right_star;
58                }
59 
60                //printf("u_left = %g u_right = %g \n",u_left, u_right);
61                soundspeed_left = sqrt(g*h_left_star);
62                soundspeed_right = sqrt(g*h_right_star);
63               
64                s_max = max(u_left+soundspeed_left, u_right+soundspeed_right);
65                if (s_max < 0.0) s_max = 0.0;
66       
67                s_min = min(u_left-soundspeed_left, u_right-soundspeed_right);
68                if (s_min > 0.0) s_min = 0.0;
69               
70               
71                // Flux formulas
72                flux_left[0] = u_left*h_left_star;
73                flux_left[1] = u_left*u_left*h_left_star + 0.5*g*h_left_star*h_left_star;
74
75                flux_right[0] = u_right*h_right_star;
76                flux_right[1] = u_right*u_right*h_right_star + 0.5*g*h_right_star*h_right_star;
77
78                // Flux computation
79                denom = s_max-s_min;
80                if (denom < epsilon && z_right > z_left) {
81                        for (i=0; i<2; i++) edgeflux[i] = 0.0;
82                        // Maximal wavespeed                   
83                        edgeflux[1] = edgeflux[1] + (0.5*g*h_left*h_left - 0.5*g*h_left_star*h_left_star );
84                        edgeflux[1] *= normals;
85                        *max_speed = 0.0;
86                } else if (denom < epsilon){
87                        for (i=0; i<2; i++) edgeflux[i] = 0.0;
88                        // Maximal wavespeed
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*(h_right_star-h_left_star);                   //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*(flux_right[0]-flux_left[0]);                 //s_max*s_min*(qr[1]-ql[1]);
96                        edgeflux[1] /= denom;
97                        edgeflux[1] = edgeflux[1] + (0.5*g*h_left*h_left - 0.5*g*h_left_star*h_left_star );
98                        edgeflux[1] *= normals;
99                        // Maximal wavespeed
100                        *max_speed = max(fabs(s_max), fabs(s_min));
101                } 
102               
103    return 0;           
104        }
105               
106               
107       
108       
109// Computational function for flux computation
110double _compute_fluxes_ext(double timestep,
111                double epsilon,
112                double g,
113               
114                long* neighbours,
115                long* neighbour_vertices,
116                double* normals,
117                double* areas,
118               
119                double* stage_edge_values,
120                double* xmom_edge_values,
121                double* bed_edge_values,
122                double* height_edge_values,
123                double* vel_edge_values,
124               
125                double* stage_boundary_values,
126                double* xmom_boundary_values,
127                double* bed_boundary_values,
128                double* height_boundary_values,
129                double* vel_boundary_values,
130               
131                double* stage_explicit_update,
132                double* xmom_explicit_update,
133                double* bed_explicit_update,
134                double* height_explicit_update,
135                double* vel_explicit_update,
136               
137                int number_of_elements,
138                double* max_speed_array) {
139               
140                double flux[2], ql[5], qr[5], edgeflux[2];
141                double max_speed, normal;
142                int k, i, ki, n, m, nm=0;
143               
144                for (k=0; k<number_of_elements; k++) {
145                        flux[0] = 0.0;
146                        flux[1] = 0.0;
147                       
148                        for (i=0; i<2; i++) {
149                                ki = k*2+i;
150                               
151                                ql[0] = stage_edge_values[ki];
152                                ql[1] = xmom_edge_values[ki];
153                                ql[2] = bed_edge_values[ki];
154                                ql[3] = height_edge_values[ki];
155                                ql[4] = vel_edge_values[ki];
156               
157                                n = neighbours[ki];
158                                printf("neighbours[ki] = %li \n",neighbours[ki]);
159                                if (n<0) {
160                                        m = -n-1;
161                                        qr[0] = stage_boundary_values[m];
162                                        qr[1] = xmom_boundary_values[m];
163                                        qr[2] = ql[2];                                                                                                                                                 
164                                        qr[3] = height_edge_values[m];                                                                                                         
165                                        qr[4] = vel_edge_values[m];                                                                                                                             
166                   
167                                } else {
168                                        m = neighbour_vertices[ki];
169                                        nm = n*2+m;
170                                        qr[0] = stage_edge_values[nm];
171                                        qr[1] = xmom_edge_values[nm];
172                                        qr[2] = bed_edge_values[nm];
173                                        qr[3] = height_edge_values[nm];                                                                                                                 
174                                        qr[4] = vel_edge_values[nm];                                                                                                                   
175                                }
176                               
177                                normal = normals[ki];
178                                _flux_function(ql, qr, normal, g, epsilon, edgeflux, &max_speed);
179                                flux[0] -= edgeflux[0];
180                                flux[1] -= edgeflux[1];
181                               
182                                // Update timestep based on edge i and possibly neighbour n
183                                if (max_speed > epsilon) {
184                        // Original CFL calculation
185                                        timestep = min(timestep, 0.5*areas[k]/max_speed); //Here, CFL=1.0 is assumed.
186                                        if (n>=0) {
187                                                timestep = min(timestep, 0.5*areas[n]/max_speed); //Here, CFL=1.0 is assumed.
188                                        }
189                                }
190            } // End edge i (and neighbour n)
191                        flux[0] /= areas[k];
192                        stage_explicit_update[k] = flux[0];
193                        flux[1] /= areas[k];
194                        xmom_explicit_update[k] = flux[1];
195                       
196                        //Keep track of maximal speeds
197                        max_speed_array[k]=max_speed;
198                }
199               
200                return timestep;       
201        }
202       
203       
204       
205
206
207
208//=========================================================================
209// Python Glue
210//=========================================================================
211PyObject *compute_fluxes_ext(PyObject *self, PyObject *args) {
212 
213    PyArrayObject
214        *neighbours, 
215        *neighbour_vertices,
216        *normals, 
217        *areas,
218       
219        *stage_edge_values,
220        *xmom_edge_values,
221        *bed_edge_values,
222        *height_edge_values,
223        *vel_edge_values,
224       
225        *stage_boundary_values,
226        *xmom_boundary_values,
227        *bed_boundary_values,
228        *height_boundary_values,
229        *vel_boundary_values,
230       
231        *stage_explicit_update,
232        *xmom_explicit_update,
233        *bed_explicit_update,
234        *height_explicit_update,
235        *vel_explicit_update,
236       
237        *max_speed_array;
238   
239  double timestep, epsilon, g;
240  int number_of_elements;
241   
242  // Convert Python arguments to C
243  if (!PyArg_ParseTuple(args, "dddOOOOOOOOOOOOOOOOOOOiO",
244                        &timestep,
245                        &epsilon,
246                        &g,
247                       
248                        &neighbours,
249                        &neighbour_vertices,
250                        &normals,
251                        &areas,
252                       
253                        &stage_edge_values,
254                        &xmom_edge_values,
255                        &bed_edge_values,
256                        &height_edge_values,
257                        &vel_edge_values,
258                       
259                        &stage_boundary_values,
260                        &xmom_boundary_values,
261                        &bed_boundary_values,
262                        &height_boundary_values,
263                        &vel_boundary_values,
264                       
265                        &stage_explicit_update,
266                        &xmom_explicit_update,
267                        &bed_explicit_update,
268                        &height_explicit_update,
269                        &vel_explicit_update,
270                       
271                        &number_of_elements,
272                        &max_speed_array)) {
273    PyErr_SetString(PyExc_RuntimeError, "swb_comp_flux_ext.c: compute_fluxes_ext could not parse input");
274    return NULL;
275  }
276
277 
278  // Call underlying flux computation routine and update
279  // the explicit update arrays
280  timestep = _compute_fluxes_ext(timestep,
281                                 epsilon,
282                                 g,
283                                 
284                                 (long*) neighbours -> data,
285                                 (long*) neighbour_vertices -> data,
286                                 (double*) normals -> data,
287                                 (double*) areas -> data,
288                                 
289                                 (double*) stage_edge_values -> data,
290                                 (double*) xmom_edge_values -> data,
291                                 (double*) bed_edge_values -> data,
292                                 (double*) height_edge_values -> data,
293                                 (double*) vel_edge_values -> data,
294                                 
295                                 (double*) stage_boundary_values -> data,
296                                 (double*) xmom_boundary_values -> data,
297                                 (double*) bed_boundary_values -> data,
298                                 (double*) height_boundary_values -> data,
299                                 (double*) vel_boundary_values -> data,
300                                 
301                                 (double*) stage_explicit_update -> data,
302                                 (double*) xmom_explicit_update -> data,
303                                 (double*) bed_explicit_update -> data,
304                                 (double*) height_explicit_update -> data,
305                                 (double*) vel_explicit_update -> data,
306                                 
307                                 number_of_elements,
308                                 (double*) max_speed_array -> data);
309
310  // Return updated flux timestep
311  return Py_BuildValue("d", timestep);
312}
313
314
315
316//-------------------------------
317// Method table for python module
318//-------------------------------
319
320static struct PyMethodDef MethodTable[] = {
321  {"compute_fluxes_ext", compute_fluxes_ext, METH_VARARGS, "Print out"},
322  {NULL, NULL}
323};
324
325// Module initialisation
326void initswb_comp_flux_ext(void){
327  Py_InitModule("swb_comp_flux_ext", MethodTable);
328  import_array();
329}
Note: See TracBrowser for help on using the repository browser.