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

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