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

Last change on this file since 7818 was 6454, checked in by steve, 15 years ago

Added Padarn's C code

File size: 14.5 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//Innermost flux function (using w=z+h)
50int _flux_function_channel(double *q_leftm,double *q_leftp, double *q_rightm,
51                           double *q_rightp, double g, double epsilon, double h0,                           double *edgeflux, double *max_speed) {
52
53    int i;
54    double flux_left[2], flux_right[2];
55    double a_leftm,w_leftm, h_leftm, d_leftm, z_leftm, u_leftm, b_leftm, soundspeed_leftm;
56     double a_leftp,w_leftp, h_leftp, d_leftp, z_leftp, u_leftp, b_leftp, soundspeed_leftp;
57 double a_rightm,w_rightm, h_rightm, d_rightm, z_rightm, u_rightm, b_rightm, soundspeed_rightm;
58 double a_rightp,w_rightp, h_rightp, d_rightp, z_rightp, u_rightp, b_rightp, soundspeed_rightp;
59    double s_maxl, s_minl,s_maxr,s_minr, denom;
60    double zphalf,zmhalf,hleftstar,hrightstar;
61    double fluxtemp1,fluxtemp0,speedtemp;
62
63    zmhalf = max(q_leftm[2],q_leftp[2]);
64    zphalf = max(q_rightm[2],q_rightp[2]);   
65
66    a_leftm  = q_leftm[0];
67    d_leftm = q_leftm[1];
68    z_leftm  = q_leftm[2];
69    h_leftm  = max(0,q_leftm[3]+q_leftm[2]-zmhalf);
70    u_leftm  = q_leftm[4];
71    b_leftm  = q_leftm[5];
72    w_leftm  = h_leftm+z_leftm;
73
74    a_leftp  = q_leftp[0];
75    d_leftp = q_leftp[1];
76    z_leftp  = q_leftp[2];
77    h_leftp  = max(0,q_leftp[3]+q_leftp[2]-zmhalf);
78    u_leftp  = q_leftp[4];
79    b_leftp  = q_leftp[5];
80    w_leftp  = h_leftp+z_leftp;
81
82    a_rightm  = q_rightm[0];
83    d_rightm = q_rightm[1];
84    z_rightm  = q_rightm[2];
85    h_rightm  = max(0,q_rightm[3]+q_rightm[2]-zphalf);
86    u_rightm  = q_rightm[4];
87    b_rightm  = q_rightm[5];
88    w_rightm  = h_rightm+z_rightm;
89
90    a_rightp  = q_rightp[0];
91    d_rightp = q_rightp[1];
92    z_rightp  = q_rightp[2];
93    h_rightp  = max(0,q_rightp[3]+q_rightp[2]-zphalf);
94    u_rightp  = q_rightp[4];
95    b_rightp  = q_rightp[5];
96    w_rightp  = h_rightp+z_rightp;
97
98    hleftstar = q_leftp[3];
99    hrightstar = q_rightm[3];   
100
101
102
103
104    soundspeed_leftp = sqrt(g*h_leftp);
105    soundspeed_leftm = sqrt(g*h_leftm);
106    soundspeed_rightp = sqrt(g*h_rightp);
107    soundspeed_rightm = sqrt(g*h_rightm);
108
109    s_maxl = max(u_leftm+soundspeed_leftm, u_leftp+soundspeed_leftp);
110    if (s_maxl < 0.0) s_maxl = 0.0;
111
112    s_minl = min(u_leftm-soundspeed_leftm, u_leftp-soundspeed_leftp);
113    if (s_minl > 0.0) s_minl = 0.0;
114
115    s_maxr = max(u_rightm+soundspeed_rightm, u_rightp+soundspeed_rightp);
116    if (s_maxr < 0.0) s_maxr = 0.0;
117
118    s_minr = min(u_rightm-soundspeed_rightm, u_rightp-soundspeed_rightp);
119    if (s_minr > 0.0) s_minr = 0.0;
120
121    // Flux formulas for left hand side
122    flux_left[0] = d_leftm;
123    flux_left[1] = u_leftm*d_leftm + 0.5*g*h_leftm*h_leftm*b_leftm;
124
125    flux_right[0] = d_leftp;
126    flux_right[1] = u_leftp*d_leftp + 0.5*g*h_leftp*h_leftp*b_leftp;
127
128   
129    // Flux computation for left hand side
130    denom = s_maxl-s_minl;
131    if (denom < epsilon) {
132        for (i=0; i<2; i++) edgeflux[i] = 0.0;
133    } else {
134        edgeflux[0] = s_maxl*flux_left[0] - s_minl*flux_right[0];
135//      edgeflux[0] += s_maxl*s_minl*(a_leftp-a_leftm);
136        edgeflux[0] /= denom;
137        edgeflux[1] = s_maxl*flux_left[1] - s_minl*flux_right[1];
138        edgeflux[1] += s_maxl*s_minl*(d_leftp-d_leftm);
139        edgeflux[1] /= denom;
140   
141
142    }
143   
144        fluxtemp0 = edgeflux[0];
145        fluxtemp1 = edgeflux[1];
146   
147
148    // Flux formulas for right hand side
149    flux_left[0] = d_rightm;
150    flux_left[1] = u_rightm*d_rightm + 0.5*g*h_rightm*h_rightm;
151
152    flux_right[0] = d_rightp;
153    flux_right[1] = u_rightp*d_rightp + 0.5*g*h_rightp*h_rightp;
154   
155   
156   
157    // Flux computation for right hand side
158    denom = s_maxr-s_minr;
159    if (denom < epsilon) {
160        for (i=0; i<2; i++) edgeflux[i] = 0.0;
161    } else {
162        edgeflux[0] = s_maxr*flux_left[0] - s_minr*flux_right[0];
163//              edgeflux[0] += s_maxr*s_minr*(a_rightp-a_rightm);
164        edgeflux[0] /= denom;
165        edgeflux[1] = s_maxr*flux_left[1] - s_minr*flux_right[1];
166        edgeflux[1] += s_maxr*s_minr*(d_rightp-d_rightm);
167        edgeflux[1] /= denom;
168   
169
170    }
171   
172   
173    edgeflux[0]=edgeflux[0]-fluxtemp0;
174    edgeflux[1]=edgeflux[1]-fluxtemp1;
175   
176   
177     
178    edgeflux[1]-=0.5*g*h_rightm*h_rightm-0.5*g*hrightstar*hrightstar+0.5*g*hleftstar*hleftstar-0.5*g*h_leftp*h_leftp;
179
180   
181
182     
183
184     //edgeflux[1]-=0.5*g*b_rightm*h_rightm*h_rightm-0.5*g*b_leftp*h_leftp*h_leftp;
185                // Maximal wavespeed
186        if ( (s_maxl-s_minl)<epsilon && (s_maxr-s_minr)<epsilon ){
187            *max_speed = 0.0;
188        }else{
189        speedtemp = max(fabs(s_maxl),fabs(s_minl));
190        speedtemp = max(speedtemp,fabs(s_maxr));
191        speedtemp = max(speedtemp,fabs(s_minr));
192        *max_speed = speedtemp;
193        }
194
195//printf("%f\n",h_right);
196    return 0;
197}
198
199
200
201
202
203// Computational function for flux computation
204double _compute_fluxes_channel_ext(double cfl,
205                               double timestep,
206                               double epsilon,
207                               double g,
208                               double h0,
209                               long* neighbours,
210                               long* neighbour_vertices,
211                               double* normals,
212                               double* areas,
213                               double* area_edge_values,
214                               double* discharge_edge_values,
215                               double* bed_edge_values,
216                               double* height_edge_values,
217                               double* velocity_edge_values,
218                               double* width_edge_values,
219                               double* area_boundary_values,
220                               double* discharge_boundary_values,
221                               double* bed_boundary_values,
222                               double* height_boundary_values,
223                               double* velocity_boundary_values,
224                               double* width_boundary_values,
225                               double* area_explicit_update,
226                               double* discharge_explicit_update,
227                               int number_of_elements,
228                               double* max_speed_array) {
229
230    double flux[2], qlm[6],qlp[6], qrm[6],qrp[6], edgeflux[2];
231    double max_speed, normal;
232    int k, i, ki, n, m, nm=0;
233    double zstar;
234    for (k=0; k<number_of_elements; k++) {
235        flux[0] = 0.0;
236        flux[1] = 0.0;
237
238       
239            ki = k*2;
240           
241
242         n = neighbours[ki];
243            if (n<0) {
244                m = -n-1;
245
246            qlm[0] = area_boundary_values[m];
247            qlm[1] = discharge_boundary_values[m];
248            qlm[2] = bed_boundary_values[m];
249            qlm[3] = height_boundary_values[m];
250            qlm[4] = velocity_boundary_values[m];
251            qlm[5] = width_boundary_values[m];
252
253            }else{
254        m = neighbour_vertices[ki];
255                nm = n*2+m;
256           
257
258            qlm[0] = area_edge_values[nm];
259            qlm[1] = discharge_edge_values[nm];
260            qlm[2] = bed_edge_values[nm];
261            qlm[3] = height_edge_values[nm];
262            qlm[4] = velocity_edge_values[nm];
263            qlm[5] = width_edge_values[nm];
264    }
265            qlp[0] = area_edge_values[ki];
266            qlp[1] = discharge_edge_values[ki];
267            qlp[2] = bed_edge_values[ki];
268            qlp[3] = height_edge_values[ki];
269            qlp[4] = velocity_edge_values[ki];
270            qlp[5] = width_edge_values[ki];
271
272         ki = k*2+1;
273           
274
275         n = neighbours[ki];
276            if (n<0) {
277                m = -n-1;
278            qrp[0] = area_boundary_values[m];
279            qrp[1] = discharge_boundary_values[m];
280            qrp[2] = bed_boundary_values[m];
281            qrp[3] = height_boundary_values[m];
282            qrp[4] = velocity_boundary_values[m];
283            qrp[5] = width_boundary_values[m];
284         
285
286
287            }else{
288        m = neighbour_vertices[ki];
289                nm = n*2+m;
290           
291           
292            qrp[0] = area_edge_values[nm];
293            qrp[1] = discharge_edge_values[nm];
294            qrp[2] = bed_edge_values[nm];
295            qrp[3] = height_edge_values[nm];
296            qrp[4] = velocity_edge_values[nm];
297            qrp[5] = width_edge_values[nm];
298            }
299            qrm[0] = area_edge_values[ki];
300            qrm[1] = discharge_edge_values[ki];
301            qrm[2] = bed_edge_values[ki];
302            qrm[3] = height_edge_values[ki];
303            qrm[4] = velocity_edge_values[ki];
304            qrm[5] = width_edge_values[ki];
305
306             _flux_function_channel(qlm,qlp,qrm,qrp,g,epsilon,h0,edgeflux,&max_speed);
307            flux[0] -= edgeflux[0];
308            flux[1] -= edgeflux[1];
309           
310            // Update timestep based on edge i and possibly neighbour n
311            if (max_speed > epsilon) {
312                // Original CFL calculation
313
314                timestep = min(timestep, 0.5*cfl*areas[k]/max_speed);
315                if (n>=0) {
316                    timestep = min(timestep, 0.5*cfl*areas[n]/max_speed);
317                }
318            }
319         // End edge i (and neighbour n)
320        flux[0] /= areas[k];
321        area_explicit_update[k] = flux[0];
322        flux[1] /= areas[k];
323        discharge_explicit_update[k] = flux[1];
324        //Keep track of maximal speeds
325        max_speed_array[k]=max_speed;
326    }
327    return timestep;
328
329}
330
331
332//-------------------------------------------------------------
333// Old code
334//------------------------------------------------------------
335//Innermost flux function (using w=z+h)
336
337
338
339
340
341// Computational function for flux computation
342
343
344//=========================================================================
345// Python Glue
346//=========================================================================
347
348
349
350//------------------------------------------------
351// New velocity based compute fluxes
352//------------------------------------------------
353
354PyObject *compute_fluxes_channel_ext(PyObject *self, PyObject *args) {
355
356    PyObject
357        *domain,
358        *area,
359        *discharge,
360        *bed,
361        *height,
362        *velocity,
363        *width;
364
365    PyArrayObject
366        *neighbours,
367        *neighbour_vertices,
368        *normals,
369        *areas,
370        *area_vertex_values,
371        *discharge_vertex_values,
372        *bed_vertex_values,
373        *height_vertex_values,
374        *velocity_vertex_values,
375        *width_vertex_values,
376        *area_boundary_values,
377        *discharge_boundary_values,
378        *bed_boundary_values,
379        *height_boundary_values,
380        *velocity_boundary_values,
381        *width_boundary_values,
382        *area_explicit_update,
383        *discharge_explicit_update,
384        *max_speed_array;
385
386  double timestep, epsilon, g, h0, cfl;
387  int number_of_elements;
388
389
390  // Convert Python arguments to C
391  if (!PyArg_ParseTuple(args, "dOOOOOOO",
392                        &timestep,
393                        &domain,
394                        &area,
395                        &discharge,
396                        &bed,
397                        &height,
398                        &velocity,
399                        &width)) {
400      PyErr_SetString(PyExc_RuntimeError, "comp_flux_channel_ext.c: compute_fluxes_channel_ext could not parse input");
401      return NULL;
402  }
403
404
405    epsilon           = get_python_double(domain,"epsilon");
406    g                 = get_python_double(domain,"g");
407    h0                = get_python_double(domain,"h0");
408    cfl               = get_python_double(domain,"CFL");
409
410
411    neighbours        = get_consecutive_array(domain, "neighbours");
412    neighbour_vertices= get_consecutive_array(domain, "neighbour_vertices");
413    normals           = get_consecutive_array(domain, "normals");
414    areas             = get_consecutive_array(domain, "areas");
415    max_speed_array   = get_consecutive_array(domain, "max_speed_array");
416
417    area_vertex_values      = get_consecutive_array(area,    "vertex_values");
418    discharge_vertex_values       = get_consecutive_array(discharge,     "vertex_values");
419    bed_vertex_values        = get_consecutive_array(bed,      "vertex_values");
420    height_vertex_values     = get_consecutive_array(height,   "vertex_values");
421    velocity_vertex_values   = get_consecutive_array(velocity, "vertex_values");
422    width_vertex_values      = get_consecutive_array(width, "vertex_values");
423
424    area_boundary_values     = get_consecutive_array(area,     "boundary_values");
425    discharge_boundary_values      = get_consecutive_array(discharge,      "boundary_values");
426    bed_boundary_values       = get_consecutive_array(bed,       "boundary_values");
427    height_boundary_values    = get_consecutive_array(height,    "boundary_values");
428    velocity_boundary_values  = get_consecutive_array(velocity,  "boundary_values");
429    width_boundary_values       = get_consecutive_array(width,     "boundary_values");
430
431
432    area_explicit_update = get_consecutive_array(area, "explicit_update");
433    discharge_explicit_update  = get_consecutive_array(discharge,  "explicit_update");
434
435    number_of_elements = area_vertex_values -> dimensions[0];
436
437    // Call underlying flux computation routine and update
438    // the explicit update arrays
439    timestep = _compute_fluxes_channel_ext(cfl,
440                                       timestep,
441                                       epsilon,
442                                       g,
443                                       h0,
444                                       (long*) neighbours -> data,
445                                       (long*) neighbour_vertices -> data,
446                                       (double*) normals -> data,
447                                       (double*) areas -> data,
448                                       (double*) area_vertex_values -> data,
449                                       (double*) discharge_vertex_values -> data,
450                                       (double*) bed_vertex_values -> data,
451                                       (double*) height_vertex_values -> data,
452                                       (double*) velocity_vertex_values -> data,
453                                       (double*) width_vertex_values -> data,
454                                       (double*) area_boundary_values -> data,
455                                       (double*) discharge_boundary_values -> data,
456                                       (double*) bed_boundary_values -> data,
457                                       (double*) height_boundary_values -> data,
458                                       (double*) velocity_boundary_values -> data,
459                                       (double*) width_boundary_values -> data,
460                                       (double*) area_explicit_update -> data,
461                                       (double*) discharge_explicit_update -> data,
462                                       number_of_elements,
463                                       (double*) max_speed_array -> data);
464
465
466    Py_DECREF(neighbours);
467    Py_DECREF(neighbour_vertices);
468    Py_DECREF(normals);
469    Py_DECREF(areas);
470    Py_DECREF(area_vertex_values);
471    Py_DECREF(discharge_vertex_values);
472    Py_DECREF(bed_vertex_values);
473    Py_DECREF(height_vertex_values);
474    Py_DECREF(velocity_vertex_values);
475    Py_DECREF(width_vertex_values);
476    Py_DECREF(area_boundary_values);
477    Py_DECREF(discharge_boundary_values);
478    Py_DECREF(bed_boundary_values);
479    Py_DECREF(height_boundary_values);
480    Py_DECREF(velocity_boundary_values);
481    Py_DECREF(width_boundary_values);
482    Py_DECREF(area_explicit_update);
483    Py_DECREF(discharge_explicit_update);
484    Py_DECREF(max_speed_array);
485
486
487    // Return updated flux timestep
488    return Py_BuildValue("d", timestep);
489}
490
491
492//-------------------------------
493// Method table for python module
494//-------------------------------
495
496static struct PyMethodDef MethodTable[] = {
497  {"compute_fluxes_channel_ext", compute_fluxes_channel_ext, METH_VARARGS, "Print out"},
498  {NULL}
499};
500
501/* // Module initialisation */
502/* void initcomp_flux_vel_ext(void){ */
503/*   Py_InitModule("comp_flux_vel_ext", MethodTable); */
504/*   import_array(); */
505/* } */
506
507void initchannel_domain_ext(void){
508  Py_InitModule("channel_domain_ext", MethodTable);
509  import_array();
510}
Note: See TracBrowser for help on using the repository browser.