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

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