source: anuga_work/development/sudi/sw_1d/periodic_waves/johns/comp_flux_ext_wellbalanced.c @ 7837

Last change on this file since 7837 was 7837, checked in by mungkasi, 12 years ago

Again, adding some codes for 1d problems on debris avalanche and periodic waves.

File size: 10.3 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_in, z_out;
18                double flux_in[2], flux_out[2];
19                double z_star, w_in, h_in, h_in_star, uh_in, soundspeed_in, u_in;
20                double w_out, h_out, h_out_star, uh_out, soundspeed_out, u_out;
21                double s_max, s_min, denom;
22               
23                w_in  = q_left[0];
24                uh_in = q_left[1];
25                uh_in = uh_in*normals;
26        z_in  = q_left[2];
27                h_in  = q_left[3];
28                u_in  = q_left[4];
29                u_in = u_in*normals;
30       
31                w_out  = q_right[0];
32                uh_out = q_right[1];
33                uh_out = uh_out*normals;
34        z_out  = q_right[2];
35                h_out  = q_right[3];
36                u_out  = q_right[4];
37                u_out = u_out*normals;
38               
39                //printf("------------------ \n");
40                //printf("w_in  = %f \n", w_in);
41                //printf("uh_in = %f \n", uh_in);
42                //printf("z_in  = %f \n", z_in);
43                //printf("h_in  = %f \n", h_in);
44                //printf("u_in  = %f \n", u_in);
45                //printf("------------------ \n");
46                //printf("w_out  = %f \n", w_out);
47                //printf("uh_out = %f \n", uh_out);
48                //printf("z_out  = %f \n", z_out);
49                //printf("h_out  = %f \n", h_out);
50                //printf("u_out  = %f \n", u_out);       
51                //////if (z_left-z_right != 0.0 ) printf("z_l - z_r = %f \n",z_left-z_right);
52               
53       
54                z_star = max(z_in, z_out);                                                                                                              //equation(7) 
55                                   
56                // Compute speeds in x-direction.
57                h_in_star = max(0.0, w_in-z_star);                                                                                                              //equation(8)
58               
59                //if (h_in_star < epsilon) {
60                //      u_in = 0.0; h_in_star = 0.0;
61                //}
62                //else {
63                //      u_in = uh_in/h_in_star;
64                //}
65                //if (h_in_star < epsilon) {
66                //h_in_star = 0.0;
67                //}
68               
69                h_out_star = max(0.0, w_out-z_star);                                                                                                    //equation(9)
70
71                //if (h_out_star < epsilon) {
72                //      u_out = 0.0; h_out_star = 0.0;
73                //}
74                //else {
75                //      u_out = uh_out/h_out_star;
76                //}
77 
78                //printf("u_left = %g u_right = %g \n",u_left, u_right);
79                soundspeed_in = sqrt(g*h_in_star);
80                soundspeed_out = sqrt(g*h_out_star);
81               
82                s_max = max(u_in+soundspeed_in, u_out+soundspeed_out);
83                if (s_max < 0.0) s_max = 0.0;
84       
85                s_min = min(u_in-soundspeed_in, u_out-soundspeed_out);
86                if (s_min > 0.0) s_min = 0.0;
87               
88               
89                // Flux formulas
90                flux_in[0] = u_in*h_in_star;
91                flux_in[1] = u_in*u_in*h_in_star + 0.5*g*h_in_star*h_in_star;
92
93                flux_out[0] = u_out*h_out_star;
94                flux_out[1] = u_out*u_out*h_out_star + 0.5*g*h_out_star*h_out_star;
95
96                // Flux computation
97                denom = s_max-s_min;
98               
99                //printf("Edge Mass Flux = %f \n", edgeflux[0]);
100                //printf("Edge Momentum Flux = %f \n", edgeflux[1]);
101               
102            if (denom < epsilon && z_out > z_in) {
103                        for (i=0; i<2; i++) edgeflux[i] = 0.0;
104                        // Maximal wavespeed           
105                        edgeflux[1] = edgeflux[1] + (0.5*g*h_in*h_in - 0.5*g*h_in_star*h_in_star );
106                        edgeflux[1] *= normals;
107                        *max_speed = 0.0;
108                        //printf("Part 1 \n");
109                } else if (denom < epsilon) {
110                        for (i=0; i<2; i++) edgeflux[i] = 0.0;
111                        // Maximal wavespeed
112                        *max_speed = 0.0;               
113                        //printf("Part 2 \n");
114                //} else if (w_out < z_in) {
115                //      for (i=0; i<2; i++) edgeflux[i] = 0.0;
116                //      *max_speed = 0.0;
117                } else {
118                        edgeflux[0] = s_max*flux_in[0] - s_min*flux_out[0];
119                        edgeflux[0] += s_max*s_min*(h_out_star-h_in_star);               
120                        edgeflux[0] /= denom;
121                        edgeflux[1] = s_max*flux_in[1] - s_min*flux_out[1];
122                        edgeflux[1] += s_max*s_min*(flux_out[0]-flux_in[0]);               
123                        edgeflux[1] /= denom;
124                        edgeflux[1] = edgeflux[1] + (0.5*g*h_in*h_in - 0.5*g*h_in_star*h_in_star);
125                        edgeflux[1] *= normals;
126                        // Maximal wavespeed
127                        *max_speed = max(fabs(s_max), fabs(s_min));
128                        //printf("Part 3 \n");
129                } 
130                //printf("Edge Mass Flux = %f \n", edgeflux[0]);
131                //printf("Edge Momentum Flux = %f \n", edgeflux[1]);
132
133    return 0;           
134        }
135               
136               
137       
138       
139// Computational function for flux computation
140double _compute_fluxes_ext_wellbalanced(double timestep,
141                double epsilon,
142                double g,
143               
144                long* neighbours,
145                long* neighbour_vertices,
146                double* normals,
147                double* areas,
148               
149                double* stage_edge_values,
150                double* xmom_edge_values,
151                double* bed_edge_values,
152                double* height_edge_values,
153                double* vel_edge_values,
154               
155                double* stage_boundary_values,
156                double* xmom_boundary_values,
157                double* bed_boundary_values,
158                double* height_boundary_values,
159                double* vel_boundary_values,
160               
161                double* stage_explicit_update,
162                double* xmom_explicit_update,
163                double* bed_explicit_update,
164                double* height_explicit_update,
165                double* vel_explicit_update,
166               
167                int number_of_elements,
168                double* max_speed_array) {
169               
170                double flux[2], ql[5], qr[5], edgeflux[2];
171                double max_speed, normal;
172                int k, i, ki, n, m, nm=0;
173               
174                for (k=0; k<number_of_elements; k++) {
175                        flux[0] = 0.0;
176                        flux[1] = 0.0;
177                       
178                        //printf("==================================== \n");
179                        //printf("==================================== \n");
180                        //printf("cell_number = %i \n",k);
181                       
182                        for (i=0; i<2; i++) {
183                               
184                                //printf("......................... \n");
185                                //printf("interface = %i \n",i);
186                               
187                                ki = k*2+i;
188                               
189                                ql[0] = stage_edge_values[ki];
190                                ql[1] = xmom_edge_values[ki];
191                                ql[2] = bed_edge_values[ki];
192                ql[3] = height_edge_values[ki];
193                                ql[4] = vel_edge_values[ki];
194               
195                                n = neighbours[ki];
196                                //printf("neighbours[ki] = %li \n",neighbours[ki]);
197                                if (n<0) {
198                                        m = -n-1;
199                                        qr[0] = stage_boundary_values[m];
200                                        qr[1] = xmom_boundary_values[m];
201                                        qr[2] = ql[2];                                                                                                                                                 
202                                        qr[3] = height_edge_values[m];                                                                                                         
203                                        qr[4] = vel_edge_values[m];                                                                                                                             
204                   
205                                } else {
206                                        m = neighbour_vertices[ki];
207                                        nm = n*2+m;
208                                        qr[0] = stage_edge_values[nm];
209                                        qr[1] = xmom_edge_values[nm];
210                                        qr[2] = bed_edge_values[nm];
211                                        qr[3] = height_edge_values[nm];                                                                                                                 
212                                        qr[4] = vel_edge_values[nm];                                                                                                                   
213                                }
214                               
215                                //Added condition!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
216                                //if (qr[0] < ql[2]) {
217                                //ql[4] = 0.0;
218                                //vel_edge_values[ki] = 0.0;
219                                //}                             
220                                //if (qr[2] > ql[2]) {
221                                //ql[1] = 0.0;
222                                //ql[4] = 0.0;
223                                //}
224                               
225                                normal = normals[ki];
226                                _flux_function_wellbalanced(ql, qr, normal, g, epsilon, edgeflux, &max_speed);
227                                flux[0] -= edgeflux[0];
228                                flux[1] -= edgeflux[1];
229                               
230                                // Update timestep based on edge i and possibly neighbour n
231                                if (max_speed > epsilon) {
232                        // Original CFL calculation
233                                        timestep = min(timestep, 0.5*areas[k]/max_speed); //Here, CFL=1.0 is assumed.
234                                        if (n>=0) {
235                                                timestep = min(timestep, 0.5*areas[n]/max_speed); //Here, CFL=1.0 is assumed.
236                                        }
237                                }
238            } // End edge i (and neighbour n)
239                        flux[0] /= areas[k];
240                        stage_explicit_update[k] = flux[0];
241                        flux[1] /= areas[k];
242                        xmom_explicit_update[k] = flux[1];
243                       
244                        //printf("......................... \n");
245                        //printf("areas = %f \n", areas[k]);
246                        //printf("Mass Flux in this cell = %f \n",flux[0]);
247                        //printf("Momentum Flux in this cell = %f \n",flux[1]);
248                        //printf("stage_explicit_update in this cell = %f \n",stage_explicit_update[k]);
249                        //printf("xmom_explicit_update in this cell = %f \n",xmom_explicit_update[k]);
250                       
251                       
252                        //Keep track of maximal speeds
253                        max_speed_array[k]=max_speed;
254                }
255               
256                return timestep;       
257        }
258       
259       
260       
261
262
263
264//=========================================================================
265// Python Glue
266//=========================================================================
267PyObject *compute_fluxes_ext_wellbalanced(PyObject *self, PyObject *args) {
268 
269    PyArrayObject
270        *neighbours, 
271        *neighbour_vertices,
272        *normals, 
273        *areas,
274       
275        *stage_edge_values,
276        *xmom_edge_values,
277        *bed_edge_values,
278        *height_edge_values,
279        *vel_edge_values,
280       
281        *stage_boundary_values,
282        *xmom_boundary_values,
283        *bed_boundary_values,
284        *height_boundary_values,
285        *vel_boundary_values,
286       
287        *stage_explicit_update,
288        *xmom_explicit_update,
289        *bed_explicit_update,
290        *height_explicit_update,
291        *vel_explicit_update,
292       
293        *max_speed_array;
294   
295  double timestep, epsilon, g;
296  int number_of_elements;
297   
298  // Convert Python arguments to C
299  if (!PyArg_ParseTuple(args, "dddOOOOOOOOOOOOOOOOOOOiO",
300                        &timestep,
301                        &epsilon,
302                        &g,
303                       
304                        &neighbours,
305                        &neighbour_vertices,
306                        &normals,
307                        &areas,
308                       
309                        &stage_edge_values,
310                        &xmom_edge_values,
311                        &bed_edge_values,
312                        &height_edge_values,
313                        &vel_edge_values,
314                       
315                        &stage_boundary_values,
316                        &xmom_boundary_values,
317                        &bed_boundary_values,
318                        &height_boundary_values,
319                        &vel_boundary_values,
320                       
321                        &stage_explicit_update,
322                        &xmom_explicit_update,
323                        &bed_explicit_update,
324                        &height_explicit_update,
325                        &vel_explicit_update,
326                       
327                        &number_of_elements,
328                        &max_speed_array)) {
329    PyErr_SetString(PyExc_RuntimeError, "comp_flux_ext_wellbalanced.c: compute_fluxes_ext_wellbalanced could not parse input");
330    return NULL;
331  }
332
333 
334  // Call underlying flux computation routine and update
335  // the explicit update arrays
336  timestep = _compute_fluxes_ext_wellbalanced(timestep,
337                                 epsilon,
338                                 g,
339                                 
340                                 (long*) neighbours -> data,
341                                 (long*) neighbour_vertices -> data,
342                                 (double*) normals -> data,
343                                 (double*) areas -> data,
344                                 
345                                 (double*) stage_edge_values -> data,
346                                 (double*) xmom_edge_values -> data,
347                                 (double*) bed_edge_values -> data,
348                                 (double*) height_edge_values -> data,
349                                 (double*) vel_edge_values -> data,
350                                 
351                                 (double*) stage_boundary_values -> data,
352                                 (double*) xmom_boundary_values -> data,
353                                 (double*) bed_boundary_values -> data,
354                                 (double*) height_boundary_values -> data,
355                                 (double*) vel_boundary_values -> data,
356                                 
357                                 (double*) stage_explicit_update -> data,
358                                 (double*) xmom_explicit_update -> data,
359                                 (double*) bed_explicit_update -> data,
360                                 (double*) height_explicit_update -> data,
361                                 (double*) vel_explicit_update -> data,
362                                 
363                                 number_of_elements,
364                                 (double*) max_speed_array -> data);
365
366  // Return updated flux timestep
367  return Py_BuildValue("d", timestep);
368}
369
370
371
372//-------------------------------
373// Method table for python module
374//-------------------------------
375
376static struct PyMethodDef MethodTable[] = {
377  {"compute_fluxes_ext_wellbalanced", compute_fluxes_ext_wellbalanced, METH_VARARGS, "Print out"},
378  {NULL, NULL}
379};
380
381// Module initialisation
382void initcomp_flux_ext_wellbalanced(void){
383  Py_InitModule("comp_flux_ext_wellbalanced", MethodTable);
384  import_array();
385}
Note: See TracBrowser for help on using the repository browser.