1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
| /** Riemann solver for St. Venant. Uses wave propagation algorithm.
Compute waves, speeds and fluctuations at each interface for the
scalar advection problem.
*/
double grav = 1;
void rp1_swe(int meqn, int mwaves, double *ql, double *qr, double *waves,
double *speeds, double *amdq, double *apdq)
{
if (meqn != 2)
{
printf("Error : meqn != 2\n");
exit(0);
}
double hl = ql[0], hr = qr[0];
double hul = ql[1], hur = qr[1];
double delta[2];
for(int m = 0; m < meqn; m++)
{
delta[m] = qr[m] - ql[m];
}
if (hl <= 0 || hr <= 0)
{
printf("hr, hl <= 0; %g %g\n", hl, hr);
exit(0);
}
double ubar, cbar, divsqrt;
divsqrt = sqrt(hl) + sqrt(hr);
ubar = (hul/sqrt(hl) + hur/sqrt(hr))/divsqrt;
cbar = sqrt(0.5*grav*(hr + hl));
double a1, a2;
a1 = 0.5*(-delta[1] + (ubar + cbar) * delta[0])/cbar;
a2 = 0.5*( delta[1] - (ubar - cbar) * delta[0])/cbar;
int mw = 0;
waves[0 + mw*meqn] = a1;
waves[1 + mw*meqn] = a1*(ubar - cbar);
speeds[mw] = ubar - cbar;
mw = 1;
waves[0 + mw*meqn] = a2;
waves[1 + mw*meqn] = a2*(ubar + cbar);
speeds[mw] = ubar + cbar;
for(int m = 0; m < meqn; m++)
{
amdq[m] = 0;
apdq[m] = 0;
}
for(int mw = 0; mw < mwaves; mw++)
for(int m = 0; m < meqn; m++)
if (speeds[mw] < 0)
amdq[m] += speeds[mw]*waves[m+mw*meqn];
else
apdq[m] += speeds[mw]*waves[m+mw*meqn];
}
|
