source: anuga_work/development/2010-projects/anuga_1d/sww/sww_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: 10.1 KB
Line 
1#include "Python.h"
2#include "numpy/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
289//-----------------------
290PyObject *compute_fluxes_ext_short(PyObject *self, PyObject *args) {
291 
292   PyObject
293        *domain,
294        *stage, 
295        *xmom, 
296        *bed;
297
298    PyArrayObject
299        *neighbours, 
300        *neighbour_vertices,
301        *normals, 
302        *areas,
303        *stage_vertex_values,
304        *xmom_vertex_values,
305        *bed_vertex_values,
306        *stage_boundary_values,
307        *xmom_boundary_values,
308        *stage_explicit_update,
309        *xmom_explicit_update,
310        *max_speed_array;
311   
312  double timestep, epsilon, g, h0, cfl;
313  int number_of_elements;
314
315   
316  // Convert Python arguments to C
317  if (!PyArg_ParseTuple(args, "dOOOO",
318                        &timestep,
319                        &domain,
320                        &stage,
321                        &xmom,
322                        &bed)) {
323      PyErr_SetString(PyExc_RuntimeError, "comp_flux_ext.c: compute_fluxes_ext_short could not parse input");
324      return NULL;
325  }
326
327
328    epsilon           = get_python_double(domain,"epsilon");
329    g                 = get_python_double(domain,"g");
330    h0                = get_python_double(domain,"h0");
331    cfl               = get_python_double(domain,"CFL");
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(
360        cfl,
361        timestep,
362        epsilon,
363        g,
364        h0,
365        (long*) neighbours -> data,
366        (long*) neighbour_vertices -> data,
367        (double*) normals -> data,
368        (double*) areas -> data,
369        (double*) stage_vertex_values -> data,
370        (double*) xmom_vertex_values -> data,
371        (double*) bed_vertex_values -> data,
372        (double*) stage_boundary_values -> data,
373        (double*) xmom_boundary_values -> data,
374        (double*) stage_explicit_update -> data,
375        (double*) xmom_explicit_update -> data,
376        number_of_elements,
377        (double*) max_speed_array -> data);
378
379
380  Py_DECREF(neighbours);
381  Py_DECREF(neighbour_vertices);
382  Py_DECREF(normals);
383  Py_DECREF(areas);
384  Py_DECREF(stage_vertex_values);
385  Py_DECREF(xmom_vertex_values);
386  Py_DECREF(bed_vertex_values);
387  Py_DECREF(stage_boundary_values);
388  Py_DECREF(xmom_boundary_values);
389  Py_DECREF(stage_explicit_update);
390  Py_DECREF(xmom_explicit_update);
391  Py_DECREF(max_speed_array);
392
393
394
395
396  // Return updated flux timestep
397  return Py_BuildValue("d", timestep);
398}
399
400 
401
402
403//-------------------------------
404// Method table for python module
405//-------------------------------
406
407static struct PyMethodDef MethodTable[] = {
408  {"compute_fluxes_ext", compute_fluxes_ext, METH_VARARGS, "Print out"},
409  {"compute_fluxes_ext_short", compute_fluxes_ext_short, METH_VARARGS, "Print out"},
410  {NULL, NULL}
411};
412
413// Module initialisation
414void initsww_comp_flux_ext(void){
415  Py_InitModule("sww_comp_flux_ext", MethodTable);
416  import_array();
417}
Note: See TracBrowser for help on using the repository browser.