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

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