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

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

Moving 2010 project

File size: 25.7 KB
Line 
1#include "Python.h"
2#include "numpy/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.