source: anuga_work/development/2010-projects/anuga_1d/channel/channel_domain_ext.c @ 7839

Last change on this file since 7839 was 7839, checked in by steve, 12 years ago

Changing name of 1d projects so that it will be easy to moveto the trunk

File size: 25.7 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//WELL BALANCED VERSION
49//Innermost flux function (using w=z+h)
50int _flux_function_channel21(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    double batemp,bphalf,bmhalf;
63
64    zmhalf = max(q_leftm[2],q_leftp[2]);
65    zphalf = max(q_rightm[2],q_rightp[2]);   
66
67    a_leftm  = q_leftm[0];
68    d_leftm = q_leftm[1];
69    z_leftm  = q_leftm[2];
70    h_leftm  = max(0,q_leftm[3]+q_leftm[2]-zmhalf);
71    u_leftm  = q_leftm[4];
72    b_leftm  = q_leftm[5];
73    w_leftm  = h_leftm+z_leftm;
74
75    a_leftp  = q_leftp[0];
76    d_leftp = q_leftp[1];
77    z_leftp  = q_leftp[2];
78    h_leftp  = max(0,q_leftp[3]+q_leftp[2]-zmhalf);
79    u_leftp  = q_leftp[4];
80    b_leftp  = q_leftp[5];
81    w_leftp  = h_leftp+z_leftp;
82
83    a_rightm  = q_rightm[0];
84    d_rightm = q_rightm[1];
85    z_rightm  = q_rightm[2];
86    h_rightm  = max(0,q_rightm[3]+q_rightm[2]-zphalf);
87    u_rightm  = q_rightm[4];
88    b_rightm  = q_rightm[5];
89    w_rightm  = h_rightm+z_rightm;
90
91    a_rightp  = q_rightp[0];
92    d_rightp = q_rightp[1];
93    z_rightp  = q_rightp[2];
94    h_rightp  = max(0,q_rightp[3]+q_rightp[2]-zphalf);
95    u_rightp  = q_rightp[4];
96    b_rightp  = q_rightp[5];
97    w_rightp  = h_rightp+z_rightp;
98
99    hleftstar = q_leftp[3];
100    hrightstar = q_rightm[3];   
101
102     bphalf = 0.5*(b_rightm+b_rightp);
103     bmhalf = 0.5*(b_leftm+b_leftp);
104     //bphalf = min(b_rightm,b_rightp);
105     //bmhalf = min(b_leftm,b_leftp);
106
107
108    soundspeed_leftp = sqrt(g*h_leftp);
109    soundspeed_leftm = sqrt(g*h_leftm);
110    soundspeed_rightp = sqrt(g*h_rightp);
111    soundspeed_rightm = sqrt(g*h_rightm);
112
113    s_maxl = max(u_leftm+soundspeed_leftm, u_leftp+soundspeed_leftp);
114    if (s_maxl < 0.0) s_maxl = 0.0;
115
116    s_minl = min(u_leftm-soundspeed_leftm, u_leftp-soundspeed_leftp);
117    if (s_minl > 0.0) s_minl = 0.0;
118
119    s_maxr = max(u_rightm+soundspeed_rightm, u_rightp+soundspeed_rightp);
120    if (s_maxr < 0.0) s_maxr = 0.0;
121
122    s_minr = min(u_rightm-soundspeed_rightm, u_rightp-soundspeed_rightp);
123    if (s_minr > 0.0) s_minr = 0.0;
124
125    // Flux formulas for left hand side
126    flux_left[0] = d_leftm;
127    flux_left[1] = u_leftm*d_leftm + 0.5*g*h_leftm*h_leftm*bmhalf;
128
129    flux_right[0] = d_leftp;
130    flux_right[1] = u_leftp*d_leftp + 0.5*g*h_leftp*h_leftp*bmhalf;
131
132   
133    // Flux computation for left hand side
134    denom = s_maxl-s_minl;
135    if (denom < epsilon) {
136        for (i=0; i<2; i++) edgeflux[i] = 0.0;
137    } else {
138        edgeflux[0] = s_maxl*flux_left[0] - s_minl*flux_right[0];
139
140        batemp = (q_leftp[3]+q_leftp[2])*bmhalf-(q_leftm[3]+q_leftm[2])*bmhalf;
141        edgeflux[0] += s_maxl*s_minl*batemp;
142        edgeflux[0] /= denom;
143        edgeflux[1] = s_maxl*flux_left[1] - s_minl*flux_right[1];
144        edgeflux[1] += s_maxl*s_minl*(d_leftp-d_leftm);
145        edgeflux[1] /= denom;
146   
147
148    }
149        fluxtemp0 = edgeflux[0];
150        fluxtemp1 = edgeflux[1];
151   
152
153    // Flux formulas for right hand side
154    flux_left[0] = d_rightm;
155    flux_left[1] = u_rightm*d_rightm + 0.5*g*h_rightm*h_rightm*bphalf;
156
157    flux_right[0] = d_rightp;
158    flux_right[1] = u_rightp*d_rightp + 0.5*g*h_rightp*h_rightp*bphalf;
159   
160   
161   
162    // Flux computation for right hand side
163    denom = s_maxr-s_minr;
164    if (denom < epsilon) {
165        for (i=0; i<2; i++) edgeflux[i] = 0.0;
166    } else {
167        edgeflux[0] = s_maxr*flux_left[0] - s_minr*flux_right[0];
168
169        batemp = (q_rightp[3]+q_rightp[2])*bphalf-(q_rightm[3]+q_rightm[2])*bphalf;
170
171        edgeflux[0] += s_maxr*s_minr*batemp;
172        edgeflux[0] /= denom;
173        edgeflux[1] = s_maxr*flux_left[1] - s_minr*flux_right[1];
174        edgeflux[1] += s_maxr*s_minr*(d_rightp-d_rightm);
175        edgeflux[1] /= denom;
176   
177
178    }
179   
180   
181    edgeflux[0]=edgeflux[0]-fluxtemp0;
182    edgeflux[1]=edgeflux[1]-fluxtemp1;
183   
184   
185    edgeflux[1]-=0.5*g*h_rightm*h_rightm*bphalf-0.5*g*hrightstar*hrightstar*b_rightm+0.5*g*hleftstar*hleftstar*b_leftp-0.5*g*h_leftp*h_leftp*bmhalf;
186
187    // printf("edgflux:%f expected:%f \n",edgeflux[1],hrightstar*hrightstar*g*0.5*b_rightm-hleftstar*hleftstar*g*0.5*b_leftp);
188
189    edgeflux[1]-=g*(1.0/6.0)*(b_rightm*(hleftstar*hleftstar+hrightstar*(hrightstar+2*z_leftp-2*z_rightm)+hleftstar*(hrightstar+z_leftp-z_rightm))-b_leftp*(hleftstar*hleftstar+hrightstar*(hrightstar-z_leftp+z_rightm)+hleftstar*(hrightstar-2*z_leftp+2*z_rightm)));
190
191   
192    //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;
193     
194     //edgeflux[1]-=0.5*g*b_rightm*h_rightm*h_rightm-0.5*g*b_leftp*h_leftp*h_leftp;
195                // Maximal wavespeed
196        if ( (s_maxl-s_minl)<epsilon && (s_maxr-s_minr)<epsilon ){
197            *max_speed = 0.0;
198        }else{
199        speedtemp = max(fabs(s_maxl),fabs(s_minl));
200        speedtemp = max(speedtemp,fabs(s_maxr));
201        speedtemp = max(speedtemp,fabs(s_minr));
202        *max_speed = speedtemp;
203        }
204
205//printf("%f\n",h_right);
206    return 0;
207}
208// GOOD BUT NOT WELL BALANCED VERSION
209int _flux_function_channel(double *q_leftm,double *q_leftp, double *q_rightm,
210                           double *q_rightp, double g, double epsilon, double h0,                           double *edgeflux, double *max_speed){
211
212    int i;
213    double flux_left[2], flux_right[2];
214    double a_leftm,w_leftm, h_leftm, d_leftm, z_leftm, u_leftm, b_leftm, soundspeed_leftm;
215     double a_leftp,w_leftp, h_leftp, d_leftp, z_leftp, u_leftp, b_leftp, soundspeed_leftp;
216 double a_rightm,w_rightm, h_rightm, d_rightm, z_rightm, u_rightm, b_rightm, soundspeed_rightm;
217 double a_rightp,w_rightp, h_rightp, d_rightp, z_rightp, u_rightp, b_rightp, soundspeed_rightp;
218    double s_maxl, s_minl,s_maxr,s_minr, denom;
219    double zphalf,zmhalf,hleftstar,hrightstar;
220    double fluxtemp1,fluxtemp0,speedtemp;
221    double batemp,bphalf,bmhalf;
222
223    zmhalf = max(q_leftm[2],q_leftp[2]);
224    zphalf = max(q_rightm[2],q_rightp[2]);   
225
226    a_leftm  = q_leftm[0];
227    d_leftm = q_leftm[1];
228    z_leftm  = q_leftm[2];
229    h_leftm  = q_leftm[3];
230    u_leftm  = q_leftm[4];
231    b_leftm  = q_leftm[5];
232    w_leftm  = h_leftm+z_leftm;
233
234    a_leftp  = q_leftp[0];
235    d_leftp = q_leftp[1];
236    z_leftp  = q_leftp[2];
237    h_leftp  = q_leftp[3];
238    u_leftp  = q_leftp[4];
239    b_leftp  = q_leftp[5];
240    w_leftp  = h_leftp+z_leftp;
241
242    a_rightm  = q_rightm[0];
243    d_rightm = q_rightm[1];
244    z_rightm  = q_rightm[2];
245    h_rightm  = q_rightm[3];
246    u_rightm  = q_rightm[4];
247    b_rightm  = q_rightm[5];
248    w_rightm  = h_rightm+z_rightm;
249
250    a_rightp  = q_rightp[0];
251    d_rightp = q_rightp[1];
252    z_rightp  = q_rightp[2];
253    h_rightp  = q_rightp[3];
254    u_rightp  = q_rightp[4];
255    b_rightp  = q_rightp[5];
256    w_rightp  = h_rightp+z_rightp;
257
258    hleftstar = q_leftp[3];
259    hrightstar = q_rightm[3];   
260
261    bphalf = min(b_rightm,b_rightp);
262    bmhalf = min(b_leftm,b_leftp);
263
264    soundspeed_leftp = sqrt(g*h_leftp);
265    soundspeed_leftm = sqrt(g*h_leftm);
266    soundspeed_rightp = sqrt(g*h_rightp);
267    soundspeed_rightm = sqrt(g*h_rightm);
268
269    s_maxl = max(u_leftm+soundspeed_leftm, u_leftp+soundspeed_leftp);
270    if (s_maxl < 0.0) s_maxl = 0.0;
271
272    s_minl = min(u_leftm-soundspeed_leftm, u_leftp-soundspeed_leftp);
273    if (s_minl > 0.0) s_minl = 0.0;
274
275    s_maxr = max(u_rightm+soundspeed_rightm, u_rightp+soundspeed_rightp);
276    if (s_maxr < 0.0) s_maxr = 0.0;
277
278    s_minr = min(u_rightm-soundspeed_rightm, u_rightp-soundspeed_rightp);
279    if (s_minr > 0.0) s_minr = 0.0;
280
281    // Flux formulas for left hand side
282    flux_left[0] = d_leftm;
283    flux_left[1] = u_leftm*d_leftm + 0.5*g*h_leftm*h_leftm*b_leftm;
284
285    flux_right[0] = d_leftp;
286    flux_right[1] = u_leftp*d_leftp + 0.5*g*h_leftp*h_leftp*b_leftp;
287
288   
289    // Flux computation for left hand side
290    denom = s_maxl-s_minl;
291    if (denom < epsilon) {
292        for (i=0; i<2; i++) edgeflux[i] = 0.0;
293    } else {
294        edgeflux[0] = s_maxl*flux_left[0] - s_minl*flux_right[0];
295
296        batemp = (q_leftp[3]+q_leftp[2])*b_leftp-(q_leftm[3]+q_leftm[2])*b_leftm;
297        edgeflux[0] += s_maxl*s_minl*batemp;
298        edgeflux[0] /= denom;
299        edgeflux[1] = s_maxl*flux_left[1] - s_minl*flux_right[1];
300        edgeflux[1] += s_maxl*s_minl*(d_leftp-d_leftm);
301        edgeflux[1] /= denom;
302   
303
304    }
305        fluxtemp0 = edgeflux[0];
306        fluxtemp1 = edgeflux[1];
307   
308
309    // Flux formulas for right hand side
310    flux_left[0] = d_rightm;
311    flux_left[1] = u_rightm*d_rightm + 0.5*g*h_rightm*h_rightm*b_rightm;
312
313    flux_right[0] = d_rightp;
314    flux_right[1] = u_rightp*d_rightp + 0.5*g*h_rightp*h_rightp*b_rightp;
315   
316   
317   
318    // Flux computation for right hand side
319    denom = s_maxr-s_minr;
320    if (denom < epsilon) {
321        for (i=0; i<2; i++) edgeflux[i] = 0.0;
322    } else {
323        edgeflux[0] = s_maxr*flux_left[0] - s_minr*flux_right[0];
324
325        batemp = (q_rightp[3]+q_rightp[2])*b_rightp-(q_rightm[3]+q_rightm[2])*b_rightm;
326
327        edgeflux[0] += s_maxr*s_minr*batemp;
328        edgeflux[0] /= denom;
329        edgeflux[1] = s_maxr*flux_left[1] - s_minr*flux_right[1];
330        edgeflux[1] += s_maxr*s_minr*(d_rightp-d_rightm);
331        edgeflux[1] /= denom;
332   
333
334    }
335   
336   
337    edgeflux[0]=edgeflux[0]-fluxtemp0;
338    edgeflux[1]=edgeflux[1]-fluxtemp1;
339   
340    edgeflux[1]-=-0.5*0.5*g*(h_rightm+h_leftp)*(b_rightm+b_leftp)*(z_rightm-z_leftp)+0.5*(h_rightm+h_leftp)*(h_rightm+h_leftp)*0.5*0.5*(b_rightm-b_leftp)*g; 
341
342    //edgeflux[1]-=0.5*g*h_rightm*h_rightm*bphalf-0.5*g*hrightstar*hrightstar*b_rightm+0.5*g*hleftstar*hleftstar*b_leftp-0.5*g*h_leftp*h_leftp*bmhalf;
343
344    // printf("edgflux:%f expected:%f \n",edgeflux[1],hrightstar*hrightstar*g*0.5*b_rightm-hleftstar*hleftstar*g*0.5*b_leftp);
345
346    //edgeflux[1]-=g*(1.0/6.0)*(b_rightm*(hleftstar*hleftstar+hrightstar*(hrightstar+2*z_leftp-2*z_rightm)+hleftstar*(hrightstar+z_leftp-z_rightm))-b_leftp*(hleftstar*hleftstar+hrightstar*(hrightstar-z_leftp+z_rightm)+hleftstar*(hrightstar-2*z_leftp+2*z_rightm)));
347
348   
349    //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;
350     
351     //edgeflux[1]-=0.5*g*b_rightm*h_rightm*h_rightm-0.5*g*b_leftp*h_leftp*h_leftp;
352                // Maximal wavespeed
353        if ( (s_maxl-s_minl)<epsilon && (s_maxr-s_minr)<epsilon ){
354            *max_speed = 0.0;
355        }else{
356        speedtemp = max(fabs(s_maxl),fabs(s_minl));
357        speedtemp = max(speedtemp,fabs(s_maxr));
358        speedtemp = max(speedtemp,fabs(s_minr));
359        *max_speed = speedtemp;
360        }
361
362//printf("%f\n",h_right);
363    return 0;
364}
365// NAIEVE VERSION
366int _flux_function_channel2(double *q_leftm,double *q_leftp, double *q_rightm,
367                           double *q_rightp, double g, double epsilon, double h0,                           double *edgeflux, double *max_speed){
368
369    int i;
370    double flux_left[2], flux_right[2];
371    double a_leftm,w_leftm, h_leftm, d_leftm, z_leftm, u_leftm, b_leftm, soundspeed_leftm;
372     double a_leftp,w_leftp, h_leftp, d_leftp, z_leftp, u_leftp, b_leftp, soundspeed_leftp;
373 double a_rightm,w_rightm, h_rightm, d_rightm, z_rightm, u_rightm, b_rightm, soundspeed_rightm;
374 double a_rightp,w_rightp, h_rightp, d_rightp, z_rightp, u_rightp, b_rightp, soundspeed_rightp;
375    double s_maxl, s_minl,s_maxr,s_minr, denom;
376    double zphalf,zmhalf,hleftstar,hrightstar;
377    double fluxtemp1,fluxtemp0,speedtemp;
378    double batemp,bphalf,bmhalf;
379
380    zmhalf = max(q_leftm[2],q_leftp[2]);
381    zphalf = max(q_rightm[2],q_rightp[2]);   
382
383    a_leftm  = q_leftm[0];
384    d_leftm = q_leftm[1];
385    z_leftm  = q_leftm[2];
386    h_leftm  = q_leftm[3];
387    u_leftm  = q_leftm[4];
388    b_leftm  = q_leftm[5];
389    w_leftm  = h_leftm+z_leftm;
390
391    a_leftp  = q_leftp[0];
392    d_leftp = q_leftp[1];
393    z_leftp  = q_leftp[2];
394    h_leftp  = q_leftp[3];
395    u_leftp  = q_leftp[4];
396    b_leftp  = q_leftp[5];
397    w_leftp  = h_leftp+z_leftp;
398
399    a_rightm  = q_rightm[0];
400    d_rightm = q_rightm[1];
401    z_rightm  = q_rightm[2];
402    h_rightm  = q_rightm[3];
403    u_rightm  = q_rightm[4];
404    b_rightm  = q_rightm[5];
405    w_rightm  = h_rightm+z_rightm;
406
407    a_rightp  = q_rightp[0];
408    d_rightp = q_rightp[1];
409    z_rightp  = q_rightp[2];
410    h_rightp  = q_rightp[3];
411    u_rightp  = q_rightp[4];
412    b_rightp  = q_rightp[5];
413    w_rightp  = h_rightp+z_rightp;
414
415    hleftstar = q_leftp[3];
416    hrightstar = q_rightm[3];   
417
418    bphalf = min(b_rightm,b_rightp);
419    bmhalf = min(b_leftm,b_leftp);
420
421    soundspeed_leftp = sqrt(g*h_leftp);
422    soundspeed_leftm = sqrt(g*h_leftm);
423    soundspeed_rightp = sqrt(g*h_rightp);
424    soundspeed_rightm = sqrt(g*h_rightm);
425
426    s_maxl = max(u_leftm+soundspeed_leftm, u_leftp+soundspeed_leftp);
427    if (s_maxl < 0.0) s_maxl = 0.0;
428
429    s_minl = min(u_leftm-soundspeed_leftm, u_leftp-soundspeed_leftp);
430    if (s_minl > 0.0) s_minl = 0.0;
431
432    s_maxr = max(u_rightm+soundspeed_rightm, u_rightp+soundspeed_rightp);
433    if (s_maxr < 0.0) s_maxr = 0.0;
434
435    s_minr = min(u_rightm-soundspeed_rightm, u_rightp-soundspeed_rightp);
436    if (s_minr > 0.0) s_minr = 0.0;
437
438    // Flux formulas for left hand side
439    flux_left[0] = d_leftm;
440    flux_left[1] = u_leftm*d_leftm + 0.5*g*h_leftm*h_leftm*b_leftm;
441
442    flux_right[0] = d_leftp;
443    flux_right[1] = u_leftp*d_leftp + 0.5*g*h_leftp*h_leftp*b_leftp;
444
445   
446    // Flux computation for left hand side
447    denom = s_maxl-s_minl;
448    if (denom < epsilon) {
449        for (i=0; i<2; i++) edgeflux[i] = 0.0;
450    } else {
451        edgeflux[0] = s_maxl*flux_left[0] - s_minl*flux_right[0];
452
453        batemp = (q_leftp[3]+q_leftp[2])*b_leftp-(q_leftm[3]+q_leftm[2])*b_leftm;
454       
455        edgeflux[0] = 0.5*(flux_left[0]+flux_right[0]);
456        edgeflux[1] = 0.5*(flux_left[1]+flux_right[1]);
457
458   
459
460    }
461        fluxtemp0 = edgeflux[0];
462        fluxtemp1 = edgeflux[1];
463   
464
465    // Flux formulas for right hand side
466    flux_left[0] = d_rightm;
467    flux_left[1] = u_rightm*d_rightm + 0.5*g*h_rightm*h_rightm*b_rightm;
468
469    flux_right[0] = d_rightp;
470    flux_right[1] = u_rightp*d_rightp + 0.5*g*h_rightp*h_rightp*b_rightp;
471   
472   
473   
474    // Flux computation for right hand side
475    denom = s_maxr-s_minr;
476    if (denom < epsilon) {
477        for (i=0; i<2; i++) edgeflux[i] = 0.0;
478    } else {
479        edgeflux[0] = s_maxr*flux_left[0] - s_minr*flux_right[0];
480
481        batemp = (q_rightp[3]+q_rightp[2])*b_rightp-(q_rightm[3]+q_rightm[2])*b_rightm;
482
483     
484        edgeflux[0] = 0.5*(flux_right[0]+flux_left[0]);
485
486
487        edgeflux[1] = 0.5*(flux_left[1]+flux_right[1]);
488   
489
490    }
491   
492   
493    edgeflux[0]=edgeflux[0]-fluxtemp0;
494    edgeflux[1]=edgeflux[1]-fluxtemp1;
495   
496    edgeflux[1]-=-0.5*0.5*g*(h_rightm+h_leftp)*(b_rightm+b_leftp)*(z_rightm-z_leftp)+0.5*(h_rightm+h_leftp)*(h_rightm+h_leftp)*0.5*0.5*(b_rightm-b_leftp)*g; 
497
498    //edgeflux[1]-=0.5*g*h_rightm*h_rightm*bphalf-0.5*g*hrightstar*hrightstar*b_rightm+0.5*g*hleftstar*hleftstar*b_leftp-0.5*g*h_leftp*h_leftp*bmhalf;
499
500    // printf("edgflux:%f expected:%f \n",edgeflux[1],hrightstar*hrightstar*g*0.5*b_rightm-hleftstar*hleftstar*g*0.5*b_leftp);
501
502    //edgeflux[1]-=g*(1.0/6.0)*(b_rightm*(hleftstar*hleftstar+hrightstar*(hrightstar+2*z_leftp-2*z_rightm)+hleftstar*(hrightstar+z_leftp-z_rightm))-b_leftp*(hleftstar*hleftstar+hrightstar*(hrightstar-z_leftp+z_rightm)+hleftstar*(hrightstar-2*z_leftp+2*z_rightm)));
503
504   
505    //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;
506     
507     //edgeflux[1]-=0.5*g*b_rightm*h_rightm*h_rightm-0.5*g*b_leftp*h_leftp*h_leftp;
508                // Maximal wavespeed
509        if ( (s_maxl-s_minl)<epsilon && (s_maxr-s_minr)<epsilon ){
510            *max_speed = 0.0;
511        }else{
512        speedtemp = max(fabs(s_maxl),fabs(s_minl));
513        speedtemp = max(speedtemp,fabs(s_maxr));
514        speedtemp = max(speedtemp,fabs(s_minr));
515        *max_speed = speedtemp;
516        }
517
518//printf("%f\n",h_right);
519    return 0;
520}
521
522
523// Computational function for flux computation
524double _compute_fluxes_channel_ext(double cfl,
525                               double timestep,
526                               double epsilon,
527                               double g,
528                               double h0,
529                               long* neighbours,
530                               long* neighbour_vertices,
531                               double* normals,
532                               double* areas,
533                               double* area_edge_values,
534                               double* discharge_edge_values,
535                               double* bed_edge_values,
536                               double* height_edge_values,
537                               double* velocity_edge_values,
538                               double* width_edge_values,
539                               double* area_boundary_values,
540                               double* discharge_boundary_values,
541                               double* bed_boundary_values,
542                               double* height_boundary_values,
543                               double* velocity_boundary_values,
544                               double* width_boundary_values,
545                               double* area_explicit_update,
546                               double* discharge_explicit_update,
547                               int number_of_elements,
548                               double* max_speed_array) {
549
550    double flux[2], qlm[6],qlp[6], qrm[6],qrp[6], edgeflux[2];
551    double max_speed, normal;
552    int k, i, ki, n, m, nm=0;
553    double zstar;
554    for (k=0; k<number_of_elements; k++) {
555        flux[0] = 0.0;
556        flux[1] = 0.0;
557
558       
559            ki = k*2;
560           
561
562         n = neighbours[ki];
563            if (n<0) {
564                m = -n-1;
565
566            qlm[0] = area_boundary_values[m];
567            qlm[1] = discharge_boundary_values[m];
568            qlm[2] = bed_boundary_values[m];
569            qlm[3] = height_boundary_values[m];
570            qlm[4] = velocity_boundary_values[m];
571            qlm[5] = width_boundary_values[m];
572
573            }else{
574        m = neighbour_vertices[ki];
575                nm = n*2+m;
576           
577
578            qlm[0] = area_edge_values[nm];
579            qlm[1] = discharge_edge_values[nm];
580            qlm[2] = bed_edge_values[nm];
581            qlm[3] = height_edge_values[nm];
582            qlm[4] = velocity_edge_values[nm];
583            qlm[5] = width_edge_values[nm];
584    }
585            qlp[0] = area_edge_values[ki];
586            qlp[1] = discharge_edge_values[ki];
587            qlp[2] = bed_edge_values[ki];
588            qlp[3] = height_edge_values[ki];
589            qlp[4] = velocity_edge_values[ki];
590            qlp[5] = width_edge_values[ki];
591
592         ki = k*2+1;
593           
594
595         n = neighbours[ki];
596            if (n<0) {
597                m = -n-1;
598            qrp[0] = area_boundary_values[m];
599            qrp[1] = discharge_boundary_values[m];
600            qrp[2] = bed_boundary_values[m];
601            qrp[3] = height_boundary_values[m];
602            qrp[4] = velocity_boundary_values[m];
603            qrp[5] = width_boundary_values[m];
604         
605
606
607            }else{
608        m = neighbour_vertices[ki];
609                nm = n*2+m;
610           
611           
612            qrp[0] = area_edge_values[nm];
613            qrp[1] = discharge_edge_values[nm];
614            qrp[2] = bed_edge_values[nm];
615            qrp[3] = height_edge_values[nm];
616            qrp[4] = velocity_edge_values[nm];
617            qrp[5] = width_edge_values[nm];
618            }
619            qrm[0] = area_edge_values[ki];
620            qrm[1] = discharge_edge_values[ki];
621            qrm[2] = bed_edge_values[ki];
622            qrm[3] = height_edge_values[ki];
623            qrm[4] = velocity_edge_values[ki];
624            qrm[5] = width_edge_values[ki];
625
626             _flux_function_channel(qlm,qlp,qrm,qrp,g,epsilon,h0,edgeflux,&max_speed);
627            flux[0] -= edgeflux[0];
628            flux[1] -= edgeflux[1];
629           
630            // Update timestep based on edge i and possibly neighbour n
631            if (max_speed > epsilon) {
632                // Original CFL calculation
633
634                timestep = min(timestep, 0.5*cfl*areas[k]/max_speed);
635                if (n>=0) {
636                    timestep = min(timestep, 0.5*cfl*areas[n]/max_speed);
637                }
638            }
639         // End edge i (and neighbour n)
640        flux[0] /= areas[k];
641        area_explicit_update[k] = flux[0];
642        flux[1] /= areas[k];
643        discharge_explicit_update[k] = flux[1];
644        //Keep track of maximal speeds
645        max_speed_array[k]=max_speed;
646    }
647    return timestep;
648
649}
650
651
652//-------------------------------------------------------------
653// Old code
654//------------------------------------------------------------
655//Innermost flux function (using w=z+h)
656
657
658
659
660
661// Computational function for flux computation
662
663
664//=========================================================================
665// Python Glue
666//=========================================================================
667
668
669
670//------------------------------------------------
671// New velocity based compute fluxes
672//------------------------------------------------
673
674PyObject *compute_fluxes_channel_ext(PyObject *self, PyObject *args) {
675
676    PyObject
677        *domain,
678        *area,
679        *discharge,
680        *bed,
681        *height,
682        *velocity,
683        *width;
684
685    PyArrayObject
686        *neighbours,
687        *neighbour_vertices,
688        *normals,
689        *areas,
690        *area_vertex_values,
691        *discharge_vertex_values,
692        *bed_vertex_values,
693        *height_vertex_values,
694        *velocity_vertex_values,
695        *width_vertex_values,
696        *area_boundary_values,
697        *discharge_boundary_values,
698        *bed_boundary_values,
699        *height_boundary_values,
700        *velocity_boundary_values,
701        *width_boundary_values,
702        *area_explicit_update,
703        *discharge_explicit_update,
704        *max_speed_array;
705
706  double timestep, epsilon, g, h0, cfl;
707  int number_of_elements;
708
709
710  // Convert Python arguments to C
711  if (!PyArg_ParseTuple(args, "dOOOOOOO",
712                        &timestep,
713                        &domain,
714                        &area,
715                        &discharge,
716                        &bed,
717                        &height,
718                        &velocity,
719                        &width)) {
720      PyErr_SetString(PyExc_RuntimeError, "comp_flux_channel_ext.c: compute_fluxes_channel_ext could not parse input");
721      return NULL;
722  }
723
724
725    epsilon           = get_python_double(domain,"epsilon");
726    g                 = get_python_double(domain,"g");
727    h0                = get_python_double(domain,"h0");
728    cfl               = get_python_double(domain,"CFL");
729
730
731    neighbours        = get_consecutive_array(domain, "neighbours");
732    neighbour_vertices= get_consecutive_array(domain, "neighbour_vertices");
733    normals           = get_consecutive_array(domain, "normals");
734    areas             = get_consecutive_array(domain, "areas");
735    max_speed_array   = get_consecutive_array(domain, "max_speed_array");
736
737    area_vertex_values      = get_consecutive_array(area,    "vertex_values");
738    discharge_vertex_values       = get_consecutive_array(discharge,     "vertex_values");
739    bed_vertex_values        = get_consecutive_array(bed,      "vertex_values");
740    height_vertex_values     = get_consecutive_array(height,   "vertex_values");
741    velocity_vertex_values   = get_consecutive_array(velocity, "vertex_values");
742    width_vertex_values      = get_consecutive_array(width, "vertex_values");
743
744    area_boundary_values     = get_consecutive_array(area,     "boundary_values");
745    discharge_boundary_values      = get_consecutive_array(discharge,      "boundary_values");
746    bed_boundary_values       = get_consecutive_array(bed,       "boundary_values");
747    height_boundary_values    = get_consecutive_array(height,    "boundary_values");
748    velocity_boundary_values  = get_consecutive_array(velocity,  "boundary_values");
749    width_boundary_values       = get_consecutive_array(width,     "boundary_values");
750
751
752    area_explicit_update = get_consecutive_array(area, "explicit_update");
753    discharge_explicit_update  = get_consecutive_array(discharge,  "explicit_update");
754
755    number_of_elements = area_vertex_values -> dimensions[0];
756
757    // Call underlying flux computation routine and update
758    // the explicit update arrays
759    timestep = _compute_fluxes_channel_ext(cfl,
760                                       timestep,
761                                       epsilon,
762                                       g,
763                                       h0,
764                                       (long*) neighbours -> data,
765                                       (long*) neighbour_vertices -> data,
766                                       (double*) normals -> data,
767                                       (double*) areas -> data,
768                                       (double*) area_vertex_values -> data,
769                                       (double*) discharge_vertex_values -> data,
770                                       (double*) bed_vertex_values -> data,
771                                       (double*) height_vertex_values -> data,
772                                       (double*) velocity_vertex_values -> data,
773                                       (double*) width_vertex_values -> data,
774                                       (double*) area_boundary_values -> data,
775                                       (double*) discharge_boundary_values -> data,
776                                       (double*) bed_boundary_values -> data,
777                                       (double*) height_boundary_values -> data,
778                                       (double*) velocity_boundary_values -> data,
779                                       (double*) width_boundary_values -> data,
780                                       (double*) area_explicit_update -> data,
781                                       (double*) discharge_explicit_update -> data,
782                                       number_of_elements,
783                                       (double*) max_speed_array -> data);
784
785
786    Py_DECREF(neighbours);
787    Py_DECREF(neighbour_vertices);
788    Py_DECREF(normals);
789    Py_DECREF(areas);
790    Py_DECREF(area_vertex_values);
791    Py_DECREF(discharge_vertex_values);
792    Py_DECREF(bed_vertex_values);
793    Py_DECREF(height_vertex_values);
794    Py_DECREF(velocity_vertex_values);
795    Py_DECREF(width_vertex_values);
796    Py_DECREF(area_boundary_values);
797    Py_DECREF(discharge_boundary_values);
798    Py_DECREF(bed_boundary_values);
799    Py_DECREF(height_boundary_values);
800    Py_DECREF(velocity_boundary_values);
801    Py_DECREF(width_boundary_values);
802    Py_DECREF(area_explicit_update);
803    Py_DECREF(discharge_explicit_update);
804    Py_DECREF(max_speed_array);
805
806
807    // Return updated flux timestep
808    return Py_BuildValue("d", timestep);
809}
810
811
812//-------------------------------
813// Method table for python module
814//-------------------------------
815
816static struct PyMethodDef MethodTable[] = {
817  {"compute_fluxes_channel_ext", compute_fluxes_channel_ext, METH_VARARGS, "Print out"},
818  {NULL}
819};
820
821/* // Module initialisation */
822/* void initcomp_flux_vel_ext(void){ */
823/*   Py_InitModule("comp_flux_vel_ext", MethodTable); */
824/*   import_array(); */
825/* } */
826
827void initchannel_domain_ext(void){
828  Py_InitModule("channel_domain_ext", MethodTable);
829  import_array();
830}
Note: See TracBrowser for help on using the repository browser.