source: anuga_work/development/anuga_1d/comp_flux_ext_wellbalanced.c @ 7818

Last change on this file since 7818 was 6694, checked in by sudi, 15 years ago

Revised codes for discontinuous bed.

File size: 8.9 KB
Line 
1#include "Python.h"
2#include "Numeric/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_wellbalanced(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_wellbalanced(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_wellbalanced(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_wellbalanced(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, "comp_flux_ext_wellbalanced.c: compute_fluxes_ext_wellbalanced 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_wellbalanced(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_wellbalanced", compute_fluxes_ext_wellbalanced, METH_VARARGS, "Print out"},
322  {NULL, NULL}
323};
324
325// Module initialisation
326void initcomp_flux_ext_wellbalanced(void){
327  Py_InitModule("comp_flux_ext_wellbalanced", MethodTable);
328  import_array();
329}
Note: See TracBrowser for help on using the repository browser.