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

Last change on this file since 5827 was 5727, checked in by steve, 16 years ago

Setting up some more unit tests for shallow_water_domain

File size: 10.1 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(double timestep,
134                           double epsilon,
135                           double g,
136                           double h0,
137                           long* neighbours,
138                           long* neighbour_vertices,
139                           double* normals,
140                           double* areas,
141                           double* stage_edge_values,
142                           double* xmom_edge_values,
143                           double* bed_edge_values,
144                           double* stage_boundary_values,
145                           double* xmom_boundary_values,
146                           double* stage_explicit_update,
147                           double* xmom_explicit_update,
148                           int number_of_elements,
149                           double* max_speed_array) {
150               
151                double flux[2], ql[2], qr[2], edgeflux[2];
152                double zl, zr, max_speed, normal;
153                int k, i, ki, n, m, nm=0;
154               
155                for (k=0; k<number_of_elements; k++) {
156                        flux[0] = 0.0;
157                        flux[1] = 0.0;
158                       
159                        for (i=0; i<2; i++) {
160                                ki = k*2+i;
161                               
162                                ql[0] = stage_edge_values[ki];
163                                ql[1] = xmom_edge_values[ki];
164                                zl = bed_edge_values[ki];
165                               
166                                n = neighbours[ki];
167                                if (n<0) {
168                                        m = -n-1;
169                                        qr[0] = stage_boundary_values[m];
170                                        qr[1] = xmom_boundary_values[m];
171                                        zr = zl;
172                                } else {
173                                        m = neighbour_vertices[ki];
174                                        nm = n*2+m;
175                                        qr[0] = stage_edge_values[nm];
176                                        qr[1] = xmom_edge_values[nm];
177                                        zr = bed_edge_values[nm];                               
178                                }
179                               
180                                normal = normals[ki];
181                                _flux_function(ql, qr, zl, zr, normal, g, epsilon, h0, edgeflux, &max_speed);
182                                flux[0] -= edgeflux[0];
183                                flux[1] -= edgeflux[1];
184                               
185                                // Update timestep based on edge i and possibly neighbour n
186                                if (max_speed > epsilon) {
187                        // Original CFL calculation
188                                       
189                                        timestep = min(timestep, 0.5*areas[k]/max_speed); //Here, CFL=1.0 is assumed. ?????????????????????????????????????????????
190                                        if (n>=0) {
191                                                timestep = min(timestep, 0.5*areas[n]/max_speed); //Here, CFL=1.0 is assumed. ?????????????????????????????????????????????
192                                        }
193                                }
194            } // End edge i (and neighbour n)
195                        flux[0] /= areas[k];
196                        stage_explicit_update[k] = flux[0];
197                        flux[1] /= areas[k];
198                        xmom_explicit_update[k] = flux[1];
199                       
200                        //Keep track of maximal speeds
201                        max_speed_array[k]=max_speed;
202                }
203                return timestep;       
204        }
205       
206       
207       
208       
209       
210
211
212
213//=========================================================================
214// Python Glue
215//=========================================================================
216PyObject *compute_fluxes_ext(PyObject *self, PyObject *args) {
217 
218    PyArrayObject
219        *neighbours, 
220        *neighbour_vertices,
221        *normals, 
222        *areas,
223        *stage_edge_values,
224        *xmom_edge_values,
225        *bed_edge_values,
226        *stage_boundary_values,
227        *xmom_boundary_values,
228        *stage_explicit_update,
229        *xmom_explicit_update,
230        *max_speed_array;
231   
232  double timestep, epsilon, g, h0;
233  int number_of_elements;
234   
235  // Convert Python arguments to C
236  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOiO",
237                        &timestep,
238                        &epsilon,
239                        &g,
240                        &h0,
241                        &neighbours,
242                        &neighbour_vertices,
243                        &normals,
244                        &areas,
245                        &stage_edge_values,
246                        &xmom_edge_values,
247                        &bed_edge_values,
248                        &stage_boundary_values,
249                        &xmom_boundary_values,
250                        &stage_explicit_update,
251                        &xmom_explicit_update,
252                        &number_of_elements,
253                        &max_speed_array)) {
254    PyErr_SetString(PyExc_RuntimeError, "comp_flux_ext.c: compute_fluxes_ext could not parse input");
255    return NULL;
256  }
257
258 
259  // Call underlying flux computation routine and update
260  // the explicit update arrays
261  timestep = _compute_fluxes_ext(timestep,
262                                 epsilon,
263                                 g,
264                                 h0,
265                                 (long*) neighbours -> data,
266                                 (long*) neighbour_vertices -> data,
267                                 (double*) normals -> data,
268                                 (double*) areas -> data,
269                                 (double*) stage_edge_values -> data,
270                                 (double*) xmom_edge_values -> data,
271                                 (double*) bed_edge_values -> data,
272                                 (double*) stage_boundary_values -> data,
273                                 (double*) xmom_boundary_values -> data,
274                                 (double*) stage_explicit_update -> data,
275                                 (double*) xmom_explicit_update -> data,
276                                 number_of_elements,
277                                 (double*) max_speed_array -> data);
278
279  // Return updated flux timestep
280  return Py_BuildValue("d", timestep);
281}
282
283
284PyObject *compute_fluxes_ext_short(PyObject *self, PyObject *args) {
285 
286   PyObject
287        *domain,
288        *stage, 
289        *xmom, 
290        *bed;
291
292    PyArrayObject
293        *neighbours, 
294        *neighbour_vertices,
295        *normals, 
296        *areas,
297        *stage_vertex_values,
298        *xmom_vertex_values,
299        *bed_vertex_values,
300        *stage_boundary_values,
301        *xmom_boundary_values,
302        *stage_explicit_update,
303        *xmom_explicit_update,
304        *max_speed_array;
305   
306  double timestep, epsilon, g, h0;
307  int number_of_elements;
308
309   
310  // Convert Python arguments to C
311  if (!PyArg_ParseTuple(args, "dOOOO",
312                        &timestep,
313                        &domain,
314                        &stage,
315                        &xmom,
316                        &bed)) {
317      PyErr_SetString(PyExc_RuntimeError, "comp_flux_ext.c: compute_fluxes_ext_short could not parse input");
318      return NULL;
319  }
320
321
322    epsilon           = get_python_double(domain,"epsilon");
323    g                 = get_python_double(domain,"g");
324    h0                = get_python_double(domain,"h0");
325 
326   
327    neighbours        = get_consecutive_array(domain, "neighbours");
328    neighbour_vertices= get_consecutive_array(domain, "neighbour_vertices"); 
329    normals           = get_consecutive_array(domain, "normals");
330    areas             = get_consecutive_array(domain, "areas");   
331    max_speed_array   = get_consecutive_array(domain, "max_speed_array");
332   
333    stage_vertex_values = get_consecutive_array(stage, "vertex_values");   
334    xmom_vertex_values  = get_consecutive_array(xmom, "vertex_values");   
335    bed_vertex_values   = get_consecutive_array(bed, "vertex_values");   
336
337    stage_boundary_values = get_consecutive_array(stage, "boundary_values");   
338    xmom_boundary_values  = get_consecutive_array(xmom, "boundary_values");   
339
340
341    stage_explicit_update = get_consecutive_array(stage, "explicit_update");   
342    xmom_explicit_update  = get_consecutive_array(xmom, "explicit_update");   
343
344
345
346    number_of_elements = stage_vertex_values -> dimensions[0];
347
348
349 
350    // Call underlying flux computation routine and update
351    // the explicit update arrays
352    timestep = _compute_fluxes_ext(timestep,
353                                   epsilon,
354                                   g,
355                                   h0,
356                                   (long*) neighbours -> data,
357                                   (long*) neighbour_vertices -> data,
358                                   (double*) normals -> data,
359                                   (double*) areas -> data,
360                                   (double*) stage_vertex_values -> data,
361                                   (double*) xmom_vertex_values -> data,
362                                   (double*) bed_vertex_values -> data,
363                                   (double*) stage_boundary_values -> data,
364                                   (double*) xmom_boundary_values -> data,
365                                   (double*) stage_explicit_update -> data,
366                                   (double*) xmom_explicit_update -> data,
367                                   number_of_elements,
368                                   (double*) max_speed_array -> data);
369
370
371  Py_DECREF(neighbours);
372  Py_DECREF(neighbour_vertices);
373  Py_DECREF(normals);
374  Py_DECREF(areas);
375  Py_DECREF(stage_vertex_values);
376  Py_DECREF(xmom_vertex_values);
377  Py_DECREF(bed_vertex_values);
378  Py_DECREF(stage_boundary_values);
379  Py_DECREF(xmom_boundary_values);
380  Py_DECREF(stage_explicit_update);
381  Py_DECREF(xmom_explicit_update);
382  Py_DECREF(max_speed_array);
383
384
385
386
387  // Return updated flux timestep
388  return Py_BuildValue("d", timestep);
389}
390
391 
392
393
394//-------------------------------
395// Method table for python module
396//-------------------------------
397
398static struct PyMethodDef MethodTable[] = {
399  {"compute_fluxes_ext", compute_fluxes_ext, METH_VARARGS, "Print out"},
400  {"compute_fluxes_ext_short", compute_fluxes_ext_short, METH_VARARGS, "Print out"},
401  {NULL, NULL}
402};
403
404// Module initialisation
405void initcomp_flux_ext(void){
406  Py_InitModule("comp_flux_ext", MethodTable);
407  import_array();
408}
Note: See TracBrowser for help on using the repository browser.