# 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; } }