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

Last change on this file since 6454 was 5832, checked in by steve, 16 years ago
File size: 8.0 KB
RevLine 
[5827]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
24
25
26//Innermost flux function (using w=z+h)
27int _flux_function_wellbalanced(double *q_left, double *q_right,
28                double normals, double g, double epsilon,
29                double *edgeflux, double *max_speed) {
30               
31                int i;
32        double z_left, z_right;
[5832]33                double flux_left[2], flux_right[2];
[5827]34                double z_star, w_left, h_left, h_left_star, uh_left, soundspeed_left, u_left;
35                double w_right, h_right, h_right_star, uh_right, soundspeed_right, u_right;                          //h_right,
36                double s_max, s_min, denom;
37               
38                w_left = q_left[0];
39                uh_left = q_left[1];
40                uh_left = uh_left*normals;
41        z_left = q_left[2];
42       
43                w_right = q_right[0];
44                uh_right = q_right[1];
45                uh_right = uh_right*normals;
46        z_right = q_right[2];
47         
48                if (z_left-z_right != 0.0 ) printf("z_l - z_r = %f \n",z_left-z_right);
49       
50                //z = (z_left+z_right)/2.0;
51                z_star = max(z_left, z_right);                                                                                                          //equation(7) 
52                                   
53                // Compute speeds in x-direction
54                h_left = w_left-z_left;                                                                         //This is used for computing the edgeflux[1].
55                h_left_star = max(0, w_left-z_star);                                                                                                            //(8)
56   
57                if (h_left_star < epsilon) {
58                        u_left = 0.0; h_left_star = 0.0;
59                }
60                else {
61                        u_left = uh_left/h_left_star;
62                }
63       
64                h_right = w_right-z_right;                                                                                                                                                                                                     
65                h_right_star = max(0, w_right-z_star);                                                                                                          //(9)
66
67                if (h_right_star < epsilon) {
68                        u_right = 0.0; h_right_star = 0.0; 
69                }
70                else {
71                        u_right = uh_right/h_right_star;
72                }
73 
74                //printf("u_left = %g u_right = %g \n",u_left, u_right);
75                soundspeed_left = sqrt(g*h_left_star);
76                soundspeed_right = sqrt(g*h_right_star);
77               
78                s_max = max(u_left+soundspeed_left, u_right+soundspeed_right);
79                if (s_max < 0.0) s_max = 0.0;
80       
81                s_min = min(u_left-soundspeed_left, u_right-soundspeed_right);
82                if (s_min > 0.0) s_min = 0.0;
83               
84               
85                // Flux formulas
86                flux_left[0] = u_left*h_left_star;
87                flux_left[1] = u_left*u_left*h_left_star + 0.5*g*h_left_star*h_left_star;
88
89                flux_right[0] = u_right*h_right_star;
90                flux_right[1] = u_right*u_right*h_right_star + 0.5*g*h_right_star*h_right_star;
91
92                // Flux computation
93                denom = s_max-s_min;
94                if (denom < epsilon) {
95                        for (i=0; i<2; i++) edgeflux[i] = 0.0;
96                        *max_speed = 0.0;
97                } else {
98                        edgeflux[0] = s_max*flux_left[0] - s_min*flux_right[0];
99                        edgeflux[0] += s_max*s_min*(h_right_star-h_left_star);                   //s_max*s_min*(qr[0]-ql[0]);
100                        edgeflux[0] /= denom;
101                        edgeflux[1] = s_max*flux_left[1] - s_min*flux_right[1];
102                        edgeflux[1] += s_max*s_min*(flux_right[0]-flux_left[0]);                 //s_max*s_min*(qr[1]-ql[1]);
103                        edgeflux[1] /= denom;
104                        edgeflux[1] *= normals;
105   
106                //edgeflux[1] = edgeflux[1] + 0.5*g*0.5*(h_left+h_right)*0.5*(h_left+h_right) - 0.5*g*h_left_star*h_left_star;
107                //edgeflux[1] = edgeflux[1] + 0.5*g*h_right*h_right - 0.5*g*h_right_star*h_right_star;
108                edgeflux[1] = edgeflux[1] + 0.5*g*h_left*h_left - 0.5*g*h_left_star*h_left_star;           //what about the h_right and h_right_star??????????????????????
109                // Maximal wavespeed
110        *max_speed = max(fabs(s_max), fabs(s_min));
111                } 
112                //printf("max_speed = %g \n",*max_speed);
113    return 0;           
114        }
115               
116               
117       
118       
119// Computational function for flux computation
120double _compute_fluxes_ext_wellbalanced(double timestep,
121                double epsilon,
122                double g,
123                long* neighbours,
124                long* neighbour_vertices,
125                double* normals,
126                double* areas,
127                double* stage_edge_values,
128                double* xmom_edge_values,
129                double* bed_edge_values,
130                double* stage_boundary_values,
131                double* xmom_boundary_values,
132                double* stage_explicit_update,
133                double* xmom_explicit_update,
134                int number_of_elements,
135                double* max_speed_array) {
136               
137                double flux[2], ql[3], qr[3], edgeflux[2];
[5832]138                double max_speed, normal;
[5827]139                int k, i, ki, n, m, nm=0;
140               
141                for (k=0; k<number_of_elements; k++) {
142                        flux[0] = 0.0;
143                        flux[1] = 0.0;
144                       
145                        for (i=0; i<2; i++) {
146                                ki = k*2+i;
147                               
148                                ql[0] = stage_edge_values[ki];
149                                ql[1] = xmom_edge_values[ki];
150                                ql[2] = bed_edge_values[ki];
151                                //ql[3] = velocity_edge_values[ki];
152                //ql[4] = height_edge_values[ki];
153               
154                                n = neighbours[ki];
155                                if (n<0) {
156                                        m = -n-1;
157                                        qr[0] = stage_boundary_values[m];
158                                        qr[1] = xmom_boundary_values[m];
159                                        qr[2] = ql[2];
160                   
161                                } else {
162                                        m = neighbour_vertices[ki];
163                                        nm = n*2+m;
164                                        qr[0] = stage_edge_values[nm];
165                                        qr[1] = xmom_edge_values[nm];
166                                        qr[2] = bed_edge_values[nm];                           
167                                }
168                               
169                                normal = normals[ki];
170                                _flux_function_wellbalanced(ql, qr, normal, g, epsilon, edgeflux, &max_speed);
171                                flux[0] -= edgeflux[0];
172                                flux[1] -= edgeflux[1];
173                               
174                                // Update timestep based on edge i and possibly neighbour n
175                                if (max_speed > epsilon) {
176                        // Original CFL calculation
177                                       
178                                        timestep = min(timestep, 0.5*areas[k]/max_speed); //Here, CFL=1.0 is assumed. ?????????????????????????????????????????????
179                                        if (n>=0) {
180                                                timestep = min(timestep, 0.5*areas[n]/max_speed); //Here, CFL=1.0 is assumed. ?????????????????????????????????????????????
181                                        }
182                                }
183            } // End edge i (and neighbour n)
184                        flux[0] /= areas[k];
185                        stage_explicit_update[k] = flux[0];
186                        flux[1] /= areas[k];
187                        xmom_explicit_update[k] = flux[1];
188                       
189                        //Keep track of maximal speeds
190                        max_speed_array[k]=max_speed;
191                }
192                //printf("timestep = %f \n",timestep);
193                return timestep;       
194        }
195       
196       
197       
198
199
200
201//=========================================================================
202// Python Glue
203//=========================================================================
204PyObject *compute_fluxes_ext_wellbalanced(PyObject *self, PyObject *args) {
205 
206    PyArrayObject
207        *neighbours, 
208        *neighbour_vertices,
209        *normals, 
210        *areas,
211        *stage_edge_values,
212        *xmom_edge_values,
213        *bed_edge_values,
214        *stage_boundary_values,
215        *xmom_boundary_values,
216        *stage_explicit_update,
217        *xmom_explicit_update,
218        *max_speed_array;
219   
220  double timestep, epsilon, g;
221  int number_of_elements;
222   
223  // Convert Python arguments to C
224  if (!PyArg_ParseTuple(args, "dddOOOOOOOOOOOiO",
225                        &timestep,
226                        &epsilon,
227                        &g,
228                        &neighbours,
229                        &neighbour_vertices,
230                        &normals,
231                        &areas,
232                        &stage_edge_values,
233                        &xmom_edge_values,
234                        &bed_edge_values,
235                        &stage_boundary_values,
236                        &xmom_boundary_values,
237                        &stage_explicit_update,
238                        &xmom_explicit_update,
239                        &number_of_elements,
240                        &max_speed_array)) {
241    PyErr_SetString(PyExc_RuntimeError, "comp_flux_ext.c: compute_fluxes_ext could not parse input");
242    return NULL;
243  }
244
245 
246  // Call underlying flux computation routine and update
247  // the explicit update arrays
248  timestep = _compute_fluxes_ext_wellbalanced(timestep,
249                                 epsilon,
250                                 g,
251                                 (long*) neighbours -> data,
252                                 (long*) neighbour_vertices -> data,
253                                 (double*) normals -> data,
254                                 (double*) areas -> data,
255                                 (double*) stage_edge_values -> data,
256                                 (double*) xmom_edge_values -> data,
257                                 (double*) bed_edge_values -> data,
258                                 (double*) stage_boundary_values -> data,
259                                 (double*) xmom_boundary_values -> data,
260                                 (double*) stage_explicit_update -> data,
261                                 (double*) xmom_explicit_update -> data,
262                                 number_of_elements,
263                                 (double*) max_speed_array -> data);
264
265  // Return updated flux timestep
266  return Py_BuildValue("d", timestep);
267}
268
269
270
271
272
273//-------------------------------
274// Method table for python module
275//-------------------------------
276
277static struct PyMethodDef MethodTable[] = {
278  {"compute_fluxes_ext_wellbalanced", compute_fluxes_ext_wellbalanced, METH_VARARGS, "Print out"},
279  {NULL, NULL}
280};
281
282// Module initialisation
283void initcomp_flux_ext_wellbalanced(void){
284  Py_InitModule("comp_flux_ext_wellbalanced", MethodTable);
285  import_array();
286}
287       
288       
Note: See TracBrowser for help on using the repository browser.