src/riemann.h

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
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
  
/* A collection of Riemann solvers for the Saint-Venant system 
   *
   * References:
   *    [1] Kurganov, A., & Levy, D. (2002). Central-upwind schemes for the
   *    Saint-Venant system. Mathematical Modelling and Numerical
   *    Analysis, 36(3), 397-425.
   */

  #define SQRT3 1.73205080756888
  #define epsilon 1e-30

  void kinetic (double hm, double hp, double um, double up, double Delta,
  	      double * fh, double * fq, double * dtmax)
  {
    double ci = sqrt(G*hm/2.);
    double Mp = max(um + ci*SQRT3, 0.);
    double Mm = max(um - ci*SQRT3, 0.);
    double cig = ci/(6.*G*SQRT3);
    *fh = cig*3.*(Mp*Mp - Mm*Mm);
    *fq = cig*2.*(Mp*Mp*Mp - Mm*Mm*Mm);
    if (Mp > 0.) {
      double dt = CFL*Delta/Mp;
      if (dt < *dtmax)
        *dtmax = dt;
    }

    ci = sqrt(G*hp/2.);
    Mp = min(up + ci*SQRT3, 0.);
    Mm = min(up - ci*SQRT3, 0.);
    cig = ci/(6.*G*SQRT3);
    *fh += cig*3.*(Mp*Mp - Mm*Mm);
    *fq += cig*2.*(Mp*Mp*Mp - Mm*Mm*Mm);
    if (Mm < - epsilon) {
      double dt = CFL*Delta/-Mm;
      if (dt < *dtmax)
        *dtmax = dt;
    }
  }

  void kurganov (double hm, double hp, double um, double up, double Delta,
  	       double * fh, double * fq, double * dtmax)
  {
    double cp = sqrt(G*hp), cm = sqrt(G*hm);
    double ap = max(up + cp, um + cm); ap = max(ap, 0.);
    double am = min(up - cp, um - cm); am = min(am, 0.);
    double qm = hm*um, qp = hp*up;
    double a = max(ap, -am);
    if (a > epsilon) {
      *fh = (ap*qm - am*qp + ap*am*(hp - hm))/(ap - am); // (4.5) of [1]
      *fq = (ap*(qm*um + G*sq(hm)/2.) - am*(qp*up + G*sq(hp)/2.) + 
  	    ap*am*(qp - qm))/(ap - am);
      double dt = CFL*Delta/a;
      if (dt < *dtmax)
        *dtmax = dt;
    }
    else
      *fh = *fq = 0.;
  }

  void hllc (double hm, double hp, double um, double up, double Delta,
  	   double * fh, double * fq, double * dtmax)
  {
    double cm = sqrt (G*hm), cp = sqrt (G*hp);
    double ustar = (um + up)/2. + cm - cp;
    double cstar = (cm + cp)/2. + (um - up)/4.;
    double SL = hm == 0. ? up - 2.*cp : min (um - cm, ustar - cstar);
    double SR = hp == 0. ? um + 2.*cm : max (up + cp, ustar + cstar);

    if (0. <= SL) {
      *fh = um*hm;
      *fq = hm*(um*um + G*hm/2.);
    }
    else if (0. >= SR) {
      *fh = up*hp;
      *fq = hp*(up*up + G*hp/2.);
    }
    else {
      double fhm = um*hm;
      double fum = hm*(um*um + G*hm/2.);
      double fhp = up*hp;
      double fup = hp*(up*up + G*hp/2.);
      *fh = (SR*fhm - SL*fhp + SL*SR*(hp - hm))/(SR - SL);
      *fq = (SR*fum - SL*fup + SL*SR*(hp*up - hm*um))/(SR - SL);
    }

    double a = max(fabs(SL), fabs(SR));
    if (a > epsilon) {
      double dt = CFL*Delta/a;
      if (dt < *dtmax)
        *dtmax = dt;
    }
  }