donat

Document Sample
donat Powered By Docstoc
					TVD Constraints and Automatic Balancing
            in Shallow Water Flows
            NumHyp, Castro-Urdiales–September 2009


                              R. Donat1

                        A. Martinez-Gavara 2

        1                 ´
            Dept. de Matematica Aplicada, U. Valencia
            2
                Dept. de E.D. y An. Numerico, U. Sevilla




                                                           Benasque09 – p.1/47
Numerical schemes for Shallow Water Flows

Motivation


     Engineering applications of Shallow Water Flows

         ocean and hydraulic engineering
         Flows in rivers, reservoirs and coastal areas

     Stationary or nearly stationary solutions need to be accurately computed.

     Numerical treatment of source terms derived from the bottom topography is a
     very delicate issue

         Source term Upwinding [Roe-86, Bermudez&Vazquez-Cendon-94]
         Well Balanced Schemes [Greenberg&Leroux,96]




                                                                           Benasque09 – p.2/47
Shallow Water Flow
Derived from the Navier-Stokes model after:
     depth averaging
     hydrostatic hypothesis
     neglecting viscosity and turbulence
     not considering wind effects nor Coriolis force.

                   System of conservation laws plus a source term:

                              Ut + F (U )x + E(U )y = S

                   0                1        0               1
      0
        h
             1       q1                               q2          0
                                                                     0
                                                                        1
               B 2                                   q1 q2
      B q1 C + B q1 + 1 gh2
      B    C                        C  B                     C   B       C
                                    C +B                       = B −ghzx C
                                    C  B                     C
      @    A   @ h
               B
                       2                              h      C   @       A
                                                   2
                    q1 q2                         q2    1
                                    A  @                     A
        q2   t                                       + gh2         −ghzy
                     h                  x         h     2     y

                                q1 = hvx ,       q2 = hvy


                                                                             Benasque09 – p.3/47
Numerical Methods - Source term upwinding

                         Ut + F (U )x + E(U )y = S(U )

 [Roe, Proc. NL. Hyp. Prob., 86]    Point-wise evaluation of S(U ) does NOT work
 [LeVeque, JCP-98]                  Fractional splitting is NOT a good idea
 [Bermudez, Vazquez, C&F, 94]       C-property
 [Greenberg, Le Roux, SINUM 96]     Well Balanced schemes




                                                                              Benasque09 – p.4/47
Numerical Methods - Source term upwinding

                                                                   Ut + F (U )x + E(U )y = S(U )

 [Roe, Proc. NL. Hyp. Prob., 86]                                                     Point-wise evaluation of S(U ) does NOT work
 [LeVeque, JCP-98]                                                                   Fractional splitting is NOT a good idea
 [Bermudez, Vazquez, C&F, 94]                                                        C-property
 [Greenberg, Le Roux, SINUM 96]                                                      Well Balanced schemes
LeVeque’s test, JCP98: Quasi-stationary flow.

                                                                                                1.0035


                                                                                                 1.003
      1.2


                                                                                                1.0025
       1                                                                                                                                fine mesh
                                                                                                                                        1st order
                                                                                                 1.002                                  2nd order
      0.8
                                                                                    Height[m]




                                                                                                                                        3rd order
                                                                                                1.0015
      0.6

                                                                                                 1.001
      0.4

                                                                                                1.0005
      0.2

                                                                                                    1
       0

                                                                                                0.9995
                                                                                                         50   100      150        200      250      300
     −0.2
            0   0.1   0.2   0.3   0.4      0.5         0.6   0.7    0.8   0.9   1                                   Position[m]
                                   bottom and initial data




                                                                                                                                                          Benasque09 – p.4/47
Numerical Methods - Source term upwinding

                          Ut + F (U )x + E(U )y = S(U )

 [Roe, Proc. NL. Hyp. Prob., 86]      Point-wise evaluation of S(U ) does NOT work
 [LeVeque, JCP-98]                    Fractional splitting is NOT a good idea
 [Bermudez, Vazquez, C&F, 94]         C-property
 [Greenberg, Le Roux, SINUM 96]       Well Balanced schemes

The spurious waves generated by schemes that are not well-balanced can distort
completely the numerical solution.
General Strategy: Combine
     conservative scheme for homogeneous conservation laws
     adequate, upwind discretization of the source term
Things to address:
     Well-Balanced, High order schemes
     Wet/Dry fronts, existing and forming.

                                                                                Benasque09 – p.4/47
Contents




           Benasque09 – p.5/47
Contents


 1. The scheme of Gascon and Corberan JCP 2001: Automatic Well Balancing




                                                                    Benasque09 – p.5/47
Contents


 1. The scheme of Gascon and Corberan JCP 2001: Automatic Well Balancing

 2. TVD-Constraints: A second order WB hybrid scheme, Martinez-Gavara, RD,
    (in preparation)




                                                                     Benasque09 – p.5/47
Contents


 1. The scheme of Gascon and Corberan JCP 2001: Automatic Well Balancing

 2. TVD-Constraints: A second order WB hybrid scheme, Martinez-Gavara, RD,
    (in preparation)

 3. The 1J-2J scheme for Balance Laws, Caselles, Haro, RD, C&F 2009




                                                                      Benasque09 – p.5/47
SC schemes for Convection Dominated Problems

                           ∂t U + ∇F (U ) = B(U , ∇U )

U = U (x, y, t)   (x, y, t) ∈ Ω×]0, T [ + appropriate initial and boundary
conditions

     Discretization on Cartesian meshes.

                                n+1   n
                               Uij − Uij
                                         + Dij = Bij
                                  ∆t




                                                                             Benasque09 – p.6/47
SC schemes for Convection Dominated Problems

                           ∂t U + ∇F (U ) = B(U , ∇U )

U = U (x, y, t)   (x, y, t) ∈ Ω×]0, T [ + appropriate initial and boundary
conditions

     Discretization on Cartesian meshes.

                                n+1   n
                               Uij − Uij
                                         + Dij = Bij
                                  ∆t


     Dij Numerical Divergence, Shock Capturing capabilities.




                                                                             Benasque09 – p.6/47
SC schemes for Convection Dominated Problems

                           ∂t U + ∇F (U ) = B(U , ∇U )

U = U (x, y, t)   (x, y, t) ∈ Ω×]0, T [ + appropriate initial and boundary
conditions

     Discretization on Cartesian meshes.

                                n+1   n
                               Uij − Uij
                                         + Dij = Bij
                                  ∆t


     Dij Numerical Divergence, Shock Capturing capabilities.
             1
                 R
     Uij ≈ |cij | c U (x, y, t)dxdy Finite Volume Schemes
                  ij




                                                                             Benasque09 – p.6/47
SC schemes for Convection Dominated Problems

                           ∂t U + ∇F (U ) = B(U , ∇U )

U = U (x, y, t)   (x, y, t) ∈ Ω×]0, T [ + appropriate initial and boundary
conditions

     Discretization on Cartesian meshes.

                                n+1   n
                               Uij − Uij
                                         + Dij = Bij
                                  ∆t


     Dij Numerical Divergence, Shock Capturing capabilities.
             1
                 R
     Uij ≈ |cij | c U (x, y, t)dxdy Finite Volume Schemes
                  ij


     Uij ≈ U (xi , yj , t). Point-Value Schemes, [Shu & Osher, JCP-86].




                                                                             Benasque09 – p.6/47
Numerical treatment of Source terms
Source Term upwinding, Roe-86                   wt + awx = s(x, w) = k(x)w (a > 0)
                               dw                            dx
Characteristic construction:      =k·w              along       =a
                               dt                            dt
                                    Z       t
         w(x, t) = w(x − at, 0) +               k(x − a(t − s))w(x − a(t − s), s)ds
                                        0


There is an upwind domain of dependence determined by a.
’Upwinding’ in the source term discretization is essential for Well Balancing (WB)




                                                                                      Benasque09 – p.7/47
Numerical treatment of Source terms
Source Term upwinding, Roe-86                   wt + awx = s(x, w) = k(x)w (a > 0)
                               dw                            dx
Characteristic construction:      =k·w              along       =a
                               dt                            dt
                                    Z       t
         w(x, t) = w(x − at, 0) +               k(x − a(t − s))w(x − a(t − s), s)ds
                                        0


There is an upwind domain of dependence determined by a.
’Upwinding’ in the source term discretization is essential for Well Balancing (WB)


Well Balanced Schemes [Greenberg & LeRoux, SINUM-96]: Schemes that
preserve steady states at the discrete level.

The C-property [Bermudez & Vazquez-Cendon, C&F-94]: A scheme satisfies the
exact C-property if it preserves exactly stationary steady states (non-moving
water). If it is not exact, but it is accurate of order O(△x2 ), it satisfies the
approximate C-property.


                                                                                      Benasque09 – p.7/47
Gascon & Corberan’s strategy

Consider the scalar equation   wt + f (w)x = s(x, w)




                                                       Benasque09 – p.8/47
Gascon & Corberan’s strategy

Consider the scalar equation                   wt + f (w)x = s(x, w)
                                                              Z   x
When w = w(x)      →    f (w)x = s(x, w)   →    f (w) = K +           s(y, w(y))dy
                                                                  ¯
                                                                  x



Construct schemes that treat the flux and any primitive of the source term
in an analogous fashion: Flux upwinding → source term upwinding




                                                                                Benasque09 – p.8/47
Gascon & Corberan’s strategy

Consider the scalar equation                                     wt + f (w)x = s(x, w)

                          Rx
Define            b=−        ¯
                            x
                                s(y, w(y, t))dy                   →      bx = −s


wt + fx = −bx     →        wt + (f + b)x = 0                 →        wt + gx = 0        g =f +b


                                                         Z   xi
                  n           n
                 gi   =   f (wi )   +   bn ;
                                         i      bn
                                                 i   =            s(y, w(y, tn ))dy,
                                                         ¯
                                                         x
                                                             Z    xi+1
             bi+1 = bi + bi,i+1 ,              bi,i+1 =                  s(y, w(y, tn ))dy
                                                                 xi




                                                                                               Benasque09 – p.8/47
Gascon & Corberan’s strategy

Consider the scalar equation                                      wt + f (w)x = s(x, w)

                           Rx
Define             b=−        ¯
                             x
                                 s(y, w(y, t))dy                   →      bx = −s


wt + fx = −bx      →        wt + (f + b)x = 0                 →        wt + gx = 0        g =f +b


                                                          Z   xi
                   n           n
                  gi   =   f (wi )   +   bn ;
                                          i      bn
                                                  i   =            s(y, w(y, tn ))dy,
                                                          ¯
                                                          x
                                                              Z    xi+1
              bi+1 = bi + bi,i+1 ,              bi,i+1 =                  s(y, w(y, tn ))dy
                                                                  xi

OBS! The computation of bi might need numerical integration,

          bi,i+1 = ˆi,i+1 + O(△xp+1 ),
      (
                   b
                                                      gi = f (wi ) + bi ≈ gi = f (wi ) + ˆi ,
                                                                          ˆ              b
          bi = ˆi + O(△xp )
               b


                                                                                                Benasque09 – p.8/47
Automatic Well Balancing

Steady (smooth) flow: ≡   wt = 0,   ⇒   fx = s = −bx   ⇔   gx = 0




                                                                   Benasque09 – p.9/47
Automatic Well Balancing

Steady (smooth) flow: ≡   wt = 0,       ⇒   fx = s = −bx   ⇔   gx = 0


                           fx      =   s   =   −bx




                                                                       Benasque09 – p.9/47
Automatic Well Balancing

Steady (smooth) flow: ≡        wt = 0,        ⇒      fx = s = −bx          ⇔   gx = 0

                     Z    xi+1          Z    xi+1        Z   xi+1
                                 fx =               s=              −bx
                         xi                 xi           xi
,




                                                                                       Benasque09 – p.9/47
Automatic Well Balancing

Steady (smooth) flow: ≡          wt = 0,          ⇒        fx = s = −bx       ⇔     gx = 0

                  Z   xi+1          Z    xi+1        Z   xi+1
    fi+1 − fi =              fx =               s=              −bx = −bi,i+1 = −(bi+1 − bi )
                  xi                    xi           xi
,




                                                                                            Benasque09 – p.9/47
Automatic Well Balancing

Steady (smooth) flow: ≡          wt = 0,          ⇒        fx = s = −bx       ⇔     gx = 0

                  Z   xi+1          Z    xi+1        Z   xi+1
    fi+1 − fi =              fx =               s=              −bx = −bi,i+1 = −(bi+1 − bi )
                  xi                    xi           xi
,
                             gi+1 = fi+1 + bi+1 = fi + bi = gi

           ut = 0       ⇒      fi+1 − fi + bi,i+1 = 0               ≡   gi+1 − gi = 0




                                                                                            Benasque09 – p.9/47
Automatic Well Balancing

Steady (smooth) flow: ≡           wt = 0,          ⇒            fx = s = −bx       ⇔    gx = 0

                   Z   xi+1          Z    xi+1        Z    xi+1
     fi+1 − fi =              fx =               s=               −bx = −bi,i+1 = −(bi+1 − bi )
                   xi                    xi               xi
,
                              gi+1 = fi+1 + bi+1 = fi + bi = gi

            ut = 0       ⇒      fi+1 − fi + bi,i+1 = 0                  ≡    gi+1 − gi = 0

The numerical flux function should be written only in terms of gi+1 − gi .

[G&C JCP-2001] Construction of a second order TVD scheme for wt + gx = 0,

                               n+1         n
                                                          ¯n            n
                                                      ˘                       ¯
                              wj     =    wj     −λ       gj+1/2   −   ¯
                                                                       gj−1/2




                                                                                                Benasque09 – p.9/47
A WB one-step, second-order scheme

         wt + gx = 0,   g(x, t) = f (w(x, t)) + b(x, t)




                                                          Benasque09 – p.10/47
A WB one-step, second-order scheme

                      wt + gx = 0,              g(x, t) = f (w(x, t)) + b(x, t)
Look for a Lax-Wendroff type second order scheme:
                                                                                           ˛n
                            n+ 1        1
                                     n+ 2                       n+ 1                     g˛
                                                                                     △t ∂ˆ ˛
    Uin+1 = Uin − λ(˜i+ 1 − gi− 1 )
                    g       ˜  2
                                                   with        ˜
                                                               gi+ 1 := gi+ 1
                                                                   2
                                                                        ˆn         +       ˛ 1.
                              2        2                           2           2     2 ∂t i+
                                                                                                  2
                                        ˛n
ˆn            n              2          ˛ 1 = (ft + bt )˛n 1 + O(△t)
                                                        ˛
gi+1/2   =   gi+1/2   + O(△x ),       ˆ
                                      gt i+              i+
                                            2                      2

                                                   Z   x                       Z   x
                                            ∂                                          ∂s
         ft = fw wt = fw (−gx ),       bt =                s(y, w(y, t))dy =              wt dy
                                            ∂t                                         ∂w




                                                                                                  Benasque09 – p.10/47
A WB one-step, second-order scheme

                      wt + gx = 0,              g(x, t) = f (w(x, t)) + b(x, t)
Look for a Lax-Wendroff type second order scheme:
                                                                                           ˛n
                            n+ 1        1
                                     n+ 2                       n+ 1                     g˛
                                                                                     △t ∂ˆ ˛
    Uin+1 = Uin − λ(˜i+ 1 − gi− 1 )
                    g       ˜  2
                                                   with        ˜
                                                               gi+ 1 := gi+ 1
                                                                   2
                                                                        ˆn         +       ˛ 1.
                              2        2                           2           2     2 ∂t i+
                                                                                                  2

                                        ˛ 1 = (ft + bt )˛n 1 + O(△t)
                                        ˛n
ˆn            n              2
                                                        ˛
gi+1/2   =   gi+1/2   + O(△x ),       ˆ
                                      gt i+              i+
                                            2                      2

                                                   Z   x                       Z   x
                                            ∂                                          ∂s
         ft = fw wt = fw (−gx ),       bt =                s(y, w(y, t))dy =              wt dy
                                            ∂t                                         ∂w

                       > g n 1 := 1 (g n + g n ),
                       8
                       > ˆi+
                       >
                       >
                       >
                       >
                              2   2 i+1        i
                       > ˛                       n
                                         ˛n gi+1 − gi   n
                             n                                  ˛n
                       <
Second order accuracy:   ˆ  ˛ 1 := −fw ˛ 1
                         gt i+            i+ 2
                                                           + bt ˛i+ 1 ,
                       >
                       >        2                  △x               2
                       >          Z x 1
                                          2 ∂s
                       > ˛             i+
                       > ˛n
                       >
                       > bt 1 =
                             i+ 2
                                                 (y, w(y, tn ))gx (y, tn )dy,
                                            ∂w
                       :
                                    0



                                                                                                  Benasque09 – p.10/47
The WB2 scheme
   Using TRAPEZOIDAL RULE
                              Z    xi+1/2
     ˛ 1 − bt ˛n 1
     ˛n       ˛
   bt i+       i−
                         =                  sw (y, w(y, tn ))gx (y, tn )dy
          2       2               xi−1/2
                              „         n    n
                                  ˛n gi+1 − gi            n    n «
                                                    ˛n gi − gi−1 △x
                         =     sw ˛i+ 1        + sw ˛i− 1            + O(△x3 )
                                      2   △x            2   △x     2
                ˛n                             ˛n
             ∂f ˛                        △t ∂s ˛
   αn 1
    i+    =λ    ˛        ,    n
                             βi+ 1     =       ˛    ,
      2      ∂w ˛ i+ 1
                                   2     2 ∂w ˛i+ 1
                     2                                     2



                       n+1    n          △t n         n        △t n
                      wi   = wi −           (Gi+ 1 − Gi− 1 ) −    Si
                                         △x      2       2     △x
                                 “                                                 ”
               n                      n          n
              Gi+ 1     =    1
                             2
                                     gi+1   +   gi   −   αn 1 (gi+1
                                                                n
                                                          i+ 2
                                                                        −    n
                                                                            gi )
                  2
                                 “                                                              ”
                 n           1        n      n            n          n      n           n
                Si      =    2
                                     βi+ 1 (gi+1     −   gi )   +   βi− 1 (gi   −      gi−1 )
                                         2                              2




                                                                                                    Benasque09 – p.11/47
Automatic Well Balancing

                                    n      n
w(x) stationary solution ⇔ wt = 0 ⇒fj+1 − fj + bn
                                                j,j+1 = 0, ∀n ≥ 0, ∀j




                                                                        Benasque09 – p.12/47
Automatic Well Balancing

                                    n      n
w(x) stationary solution ⇔ wt = 0 ⇒fj+1 − fj + bn
                                                j,j+1 = 0, ∀n ≥ 0, ∀j

                                    R xi+1
If    n    n
     gj = fj + bn ( ≡ bi,i+1 =
                j                    xi
                                              s)       ⇒

                      n      n    n      n
                     gj+1 − gj = fj+1 − fj + bi,i+1 = 0, ∀j, ∀n

                              “                                               ”
              n                    n          n
             Gi+ 1   =    1
                          2
                                  gi+1   +   gi   −   αn 1 (gi+1
                                                             n
                                                       i+ 2
                                                                   −
                                                                 = gin  n
                                                                       gi )
                 2
                              “                                        ”
                n         1      n      n      n      n      n    n
               Si    =    2
                                βi+ 1 (gi+1 − gi ) + βi− 1 (gi − gi−1 ) = 0
                                     2                             2




                                                                                  Benasque09 – p.12/47
Automatic Well Balancing

                                    n      n
w(x) stationary solution ⇔ wt = 0 ⇒fj+1 − fj + bn
                                                j,j+1 = 0, ∀n ≥ 0, ∀j

                                    R xi+1
If    n    n
     gj = fj + bn ( ≡ bi,i+1 =
                j                    xi
                                              s)       ⇒

                      n      n    n      n
                     gj+1 − gj = fj+1 − fj + bi,i+1 = 0, ∀j, ∀n

                              “                                               ”
              n                    n          n
             Gi+ 1   =    1
                          2
                                  gi+1   +   gi   −   αn 1 (gi+1
                                                             n
                                                       i+ 2
                                                                   −
                                                                 = gin  n
                                                                       gi )
                 2
                              “                                        ”
                n         1      n      n      n      n      n    n
               Si    =    2
                                βi+ 1 (gi+1 − gi ) + βi− 1 (gi − gi−1 ) = 0
                                     2                             2



          n+1    n       △t n         n        △t n
         wi   = wi −        (Gi+ 1 − Gi− 1 ) −    Si = Uin                        exactly WB
                         △x      2       2     △x




                                                                                               Benasque09 – p.12/47
Automatic Well Balancing

                                    n      n
w(x) stationary solution ⇔ wt = 0 ⇒fj+1 − fj + bn
                                                j,j+1 = 0, ∀n ≥ 0, ∀j


OBS! When the computation of bj,j+1 needs numerical integration,

         bj,j+1 = ˆj,j+1 + O(△xp+1 ),
     (
                  b
                                         gj = f (wj ) + bj ≈ gj = f (wj ) + ˆj ,
                                                             ˆ              b
         bj = ˆj + O(△x )
              b         p



           gj+1 − gj
           ˆ      ˆ    =    fj+1 − fj + ˆj,j+1 = gj+1 − gj + O(△xp+1 )
                                        b




                                                                               Benasque09 – p.13/47
Automatic Well Balancing

                                    n      n
w(x) stationary solution ⇔ wt = 0 ⇒fj+1 − fj + bn
                                                j,j+1 = 0, ∀n ≥ 0, ∀j


OBS! When the computation of bj,j+1 needs numerical integration,

           bj,j+1 = ˆj,j+1 + O(△xp+1 ),
     (
                    b
                                              gj = f (wj ) + bj ≈ gj = f (wj ) + ˆj ,
                                                                  ˆ              b
           bj = ˆj + O(△x )
                b            p



             gj+1 − gj
             ˆ      ˆ     =      fj+1 − fj + ˆj,j+1 = gj+1 − gj + O(△xp+1 )
                                             b

                     “                                ”
    n
   Gi+ 1     =   1
                 2
                        n        n
                                                             gn   ˆn
                      gi+1 + gi − αi+ 1 (ˆi+1 − gi ) = 1 (ˆi+1 + gi ) + O(△xp+1 )
                       ˆ       ˆ     n
                                          g n
                                                  ˆ n
                                                           2
       2                                2
                     “                                       ”
      n
     Si      =   1
                 2
                      βi+ 1 (ˆi+1 − gi ) + βi− 1 (ˆi − gi−1 ) = O(△xp+1 )
                         n
                             g n
                                    ˆn       n
                                                  g n
                                                        ˆn
                         2                    2




                                                                                    Benasque09 – p.13/47
Automatic Well Balancing

                                    n      n
w(x) stationary solution ⇔ wt = 0 ⇒fj+1 − fj + bn
                                                j,j+1 = 0, ∀n ≥ 0, ∀j


OBS! When the computation of bj,j+1 needs numerical integration,

             bj,j+1 = ˆj,j+1 + O(△xp+1 ),
        (
                      b
                                                 gj = f (wj ) + bj ≈ gj = f (wj ) + ˆj ,
                                                                     ˆ              b
             bj = ˆj + O(△x )
                  b             p



               gj+1 − gj
               ˆ      ˆ      =      fj+1 − fj + ˆj,j+1 = gj+1 − gj + O(△xp+1 )
                                                b

                        “                                ”
    n
   Gi+ 1       =    1
                    2
                           n        n
                                                                gn   ˆn
                         gi+1 + gi − αi+ 1 (ˆi+1 − gi ) = 1 (ˆi+1 + gi ) + O(△xp+1 )
                          ˆ       ˆ     n
                                             g n
                                                     ˆ n
                                                              2
       2                                   2
                        “                                       ”
         n
        Si     =    1
                    2
                         βi+ 1 (ˆi+1 − gi ) + βi− 1 (ˆi − gi−1 ) = O(△xp+1 )
                            n
                                g n
                                       ˆn       n
                                                     g n
                                                           ˆn
                            2                    2


 n       n
Gi+ 1 − Gi− 1 = O(△xp+1 )                n
                                        Si = O(△xp+1 )      ⇒     n+1
                                                                 wi   = wi + O(△xp+1 )
                                                                         n
    2           2


The scheme satisfies the approximate C-property (p ≥ 2)


                                                                                       Benasque09 – p.13/47
FSB and Automatic WB


Incorporating the source term into the flux divergence ⇒

Flux-Gradient and Source-term Balancing FSB ⇒ Automatic WB

                                                          R xi+1
FSB : wt = 0 → fi+1 − fi + bi,i+1 = 0,    ∀i , bi,i+1 =    xi
                                                                   s. Key to WB

                    △x
In general ˆi,i+1 =
           b           (si + si+1 ) + O(△x3 )
                     2
Trapezoidal rule ⇒
     Second order accuracy OK
     wt = 0 → fi+1 − fi + ˆi,i+1 = O(△x3 ),
                          b                     ∀i Approximate C-property.




                                                                                  Benasque09 – p.14/47
Embid’s Problem
Scalar model for gas flow through a duct of variable cross section.
                  (
                      wt + w2 /2 x = (6x − 3)w,
                            `     ´
                                                                                    0<x<1
                      w(0, t) = 1; w(1, t) = −.1

                                                  (
                                                        1 + 3x2 − 3x,                         x < .18
Initial data: Steady solution w(x) =
                                                        −0.1 + 3x2 − 3x,                      x > .18

                           1


                          0.8


                          0.6


                          0.4


                          0.2


                           0


                         −0.2


                         −0.4


                         −0.6


                         −0.8


                          −1
                                0   0.1   0.2   0.3   0.4   0.5   0.6   0.7   0.8   0.9   1



                       WB2-scheme. Numerical steady state.
                                                                                                        Benasque09 – p.15/47
Greenberg&Leroux tests [SINUM,96]
                                                                             (
                          2
                         u                                                       1.3,         x<0
                  ut +     = −ax (x)u,                u0 (x) + a(x) =
                         2                                                       1            x>0

                                              the WB2 scheme

 1.4                               1.4                                1.5




 1.2                               1.2




  1                                 1
                                                                       1


 0.8                               0.8




 0.6                               0.6


                                                                      0.5
 0.4                               0.4




 0.2                               0.2




  0                                 0                                  0
 −0.5   0   0.5     1    1.5   2   −0.5   0     0.5     1   1.5   2   −0.5       0      0.5    1    1.5   2




                           a(x), numerical solution, hx = .02, λ = .5




                                                                                                              Benasque09 – p.16/47
Computing bi,i+1 , ˆi,i+1 in Embid’s, Greenberg-Leroux tests.
                   b

                                          u2
                                    ut + ( )x = −ax (x)u
                                          2
                  Z   xi+1
                                            △x `
                                                 ax u|xi+1 + ax u|xi + O(△x3 )
                                                                    ´
       bi,i+1 =              ax (x)u(x)dx =
                   xi                        2

                                        „                                        «
                                        ai+1 − ai    △x
           ax u|xi+1         =   ui+1             −     (uxx )i + O(△x2 )
                                           △x         2
                                    „                                   «
                                      ai+1 − ai   △x
             ax u|xi         =   ui             +     (uxx )i + O(△x2 )
                                         △x        2
                                    „               «
                               △x       ai+1 − ai
                  bi,i+1     =                          (ui+1 + ui ) + O(△x3 )
                                2          △x

                              ˆi,i+1 = 1 (ai+1 − ai ) (ui+1 + ui )
                              b
                                       2



                                                                                     Benasque09 – p.17/47
Computing bi,i+1 , ˆi,i+1 in Embid’s, Greenberg-Leroux tests.
                   b

                                    u2
                             ut + ( )x = −ax (x)u
                                     2
            (                            (
              Steady state solutions       uj+1 − uj = aj+1 − aj
                                       →
              u+a=C                        ux = −ax


                       xi+1                         xi+1                         xi+1
                                                                                         u2
                   Z                           Z                             Z
   bi,i+1     =               ax (x)u(x)dx =               −ux (x)u(x)dx =              − dx
                     xi                            xi                        xi          2
                    1                          1
              =    − (ui+1 − ui )(ui+1 + ui ) = (ai+1 − ai )(ui+1 + ui ) = ˆi,i+1
                                                                           b
                    2                          2

            gj+1 − gj = fi+1 − fi + ˆi,i+1 = fi+1 − fi + bi,i+1 = gj+1 − gj
            ˆ      ˆ                b


       ⇒          The WB2 scheme is Exactly WB for Steady States u + a = C




                                                                                           Benasque09 – p.18/47
Greenberg&Leroux tests [SINUM,96]
                                                                        (
                                  2
                                 u                                                1.3,       x<0
                         ut +      = −ax (x)u,               u+a=
                                 2                                                1          x>0

                                                 the WB2 scheme
 1.4                                  1.4                                   1.5




 1.2                                  1.2




  1                                    1
                                                                             1


 0.8                                  0.8




 0.6                                  0.6


                                                                            0.5
 0.4                                  0.4




 0.2                                  0.2




  0                                    0                                     0
 −0.5   0   0.5      1     1.5    2   −0.5   0     0.5   1    1.5   2       −0.5         0   0.5   1   1.5   2


 1.4                                  1.4                                   1.4




 1.2                                  1.2                                   1.2




  1                                    1                                     1




 0.8                                  0.8                                   0.8




 0.6                                  0.6                                   0.6




 0.4                                  0.4                                   0.4




 0.2                                  0.2                                   0.2




  0                                    0                                     0
 −0.5   0   0.5      1     1.5    2   −0.5   0     0.5   1    1.5   2       −0.5         0   0.5   1   1.5   2




                  Curbing the oscillations: A WB second order hybrid scheme

                                                                                                                 Benasque09 – p.19/47
Curbing the oscillations: G & C JCP-01
                                             n+1        n
                                                                     ˆn            n
                                                                 `                       ´
             wt + gx = 0               →    wj     =   wj   −λ       gj+1/2   −   ˆ
                                                                                  gj−1/2

G&C start with Lax-Wendroff type flux
                                                    ˛n                  !
                           1                     ∂g ˛
                ˆn
                gj+1/2 =            n    n
                                   gj + gj+1 − λ    ˛        n      n
                                                           (gj+1 − gj )
                           2                     ∂w ˛j+1/2


                               ∂g   ∂f   ∂b                              ∂f              ∂b
           g = f + b,             =    +    ,                α≈λ            ,      β≈
                               ∂w   ∂w   ∂w                              ∂w              ∂w
                8   fj+1 −fj
                                                        8
                                                            bj+1 −bj
                <   wj+1 −wj
                               ,                        <
                                                            wj+1 −wj
                                                                         ,        if wj+1 − wj = 0
   αj+1/2 = λ                              βj+1/2 = λ
                    ∂f ˛
                       ˛
                :
                    ∂w j+1/2
                                   ,                    : 0,                      if wj+1 − wj = 0




                                                                                                Benasque09 – p.20/47
Curbing the oscillations: G&C TVD scheme

               8   fj+1 −fj
                                                   8
                                                       bj+1 −bj
               <   wj+1 −wj
                              ,                    <
                                                       wj+1 −wj
                                                                  ,   if wj+1 − wj = 0
  αj+1/2 = λ                          βj+1/2 = λ
                   ∂f ˛
                      ˛
               :
                   ∂w j+1/2
                                  ,                : 0,               if wj+1 − wj = 0




                                                                                    Benasque09 – p.21/47
Curbing the oscillations: G&C TVD scheme

                8   fj+1 −fj
                                                    8
                                                        bj+1 −bj
                <   wj+1 −wj
                               ,                    <
                                                        wj+1 −wj
                                                                   ,   if wj+1 − wj = 0
   αj+1/2 = λ                          βj+1/2 = λ
                    ∂f ˛
                       ˛
                :
                    ∂w j+1/2
                                   ,                : 0,               if wj+1 − wj = 0

Use Harten’s (JCP-97) construction of TVD schemes:

                       1`                                             ´
            ¯
            gj+1/2   =    gj + gj+1 − h(αj+1/2 + βj+1/2 )(gj+1 − gj )
                       2

                            n+1    n
Adjust h(x) so that scheme wj   = wj − λ(¯j+1/2 − gj−1/2 )
                                         g        ¯                         is TVD


→ CFL-like restrictions on αi+1/2 + βi+1/2

First order scheme!. Upgrade to second order following Harten JCP-97: Apply the
first order TVD scheme to a modified PDE.




                                                                                     Benasque09 – p.21/47
Curbing the oscillations: G&C TVD scheme

                8   fj+1 −fj
                                                    8
                                                        bj+1 −bj
                <   wj+1 −wj
                               ,                    <
                                                        wj+1 −wj
                                                                   ,   if wj+1 − wj = 0
   αj+1/2 = λ                          βj+1/2 = λ
                    ∂f ˛
                       ˛
                :
                    ∂w j+1/2
                                   ,                : 0,               if wj+1 − wj = 0

Use Harten’s (JCP-97) construction of TVD schemes:

                       1`                                             ´
            ¯
            gj+1/2   =    gj + gj+1 − h(αj+1/2 + βj+1/2 )(gj+1 − gj )
                       2

                            n+1    n
Adjust h(x) so that scheme wj   = wj − λ(¯j+1/2 − gj−1/2 )
                                         g        ¯                         is TVD


→ CFL-like restrictions on αi+1/2 + βi+1/2

First order scheme!. Upgrade to second order following Harten JCP-97: Apply the
first order TVD scheme to a modified PDE.

Applications: Nozze-flows


                                                                                     Benasque09 – p.21/47
Revisiting Gascon and Corberan TVD scheme
Anna Martinez Gavara PhD Thesis.

BUT!!

     Solutions to scalar balance laws might not be TVD! TVD restrictions might be
     too stringent on numerical schemes for general balance laws.

     CFL-like restrictions based on αi+1/2 + βi+1/2 might be too restrictive for
     general balance laws




                                                                             Benasque09 – p.22/47
A Well Balanced second order Hybrid scheme
                                             △t                     △t n
Second order method          Uin+1 = Uin −   △x
                                                  n
                                                (Gi+ 1     n
                                                       − Gi− 1 ) − △x Si
                                                     2       2
                                                                            ˛n
                         “                              ”                ∂f ˛
    n       HI
   Gi+ 1 = Gi+ 1 =   1
                     2
                            n      n
                                        i+ 2
                                              n      n
                           gi+1 + gi − αn 1 (gi+1 − gi ) ,      αn 1 = λ ˛
                                                                  i+ 2
       2      2                                                          ∂u ˛i+ 1
                                                                                2




                                                                               Benasque09 – p.23/47
A Well Balanced second order Hybrid scheme
                                              △t                     △t n
Second order method           Uin+1 = Uin −   △x
                                                   n
                                                 (Gi+ 1     n
                                                        − Gi− 1 ) − △x Si
                                                      2       2
                                                                             ˛n
                          “                              ”                ∂f ˛
     n       HI
    Gi+ 1 = Gi+ 1 =   1
                      2
                             n      n
                                         i+ 2
                                               n      n
                            gi+1 + gi − αn 1 (gi+1 − gi ) ,      αn 1 = λ ˛
                                                                   i+ 2
        2       2                                                         ∂u ˛i+ 1
                                                                                 2



Lax-Wendroff-type scheme: Oscillatory behavior at discontinuous fronts

Apply the Flux-Limiter technology to curb oscillations.




                                                                                Benasque09 – p.23/47
A Well Balanced second order Hybrid scheme
                                                        △t                  △t n
Second order method                  Uin+1 = Uin −      △x
                                                             n
                                                           (Gi+ 1  n
                                                               − Gi− 1 ) − △x Si
                                                                2    2
                                                                                    ˛n
                                 “                              ”                ∂f ˛
    n       HI
   Gi+ 1 = Gi+ 1 =           1
                             2
                                    n      n
                                                i+ 2
                                                      n      n
                                   gi+1 + gi − αn 1 (gi+1 − gi ) ,      αn 1 = λ ˛
                                                                          i+ 2
       2             2                                                           ∂u ˛i+ 1
                                                                                         2



                                  n       LO             HI      LO
                                 Gi+ 1 = Gi+ 1 + φi+ 1 (Gi+ 1 − Gi+ 1 )
                                     2         2        2      2     2

                    “                                     ”
     LO         1     n      n          n      n      n
    Gi+ 1   =   2
                     gi+1 + gi − sign(αi+ 1 )(gi+1 − gi )
        2                                  2
                    “                               ”
     HI         1     n      n     n     n      n
    Gi+ 1 =     2
                     gi+1 + gi − αi+ 1 (gi+1 − gi )
        2                                  2

    φi+ 1 = φ(ri+ 1 )
        2                2

                                                         > gi − gi−1 ,
                                                         8
                                                         >
                                                         < g               αi+ 1 > 0;
                                                             i+1 − gi         2
    Upwind smoothness indicator ri+ 1                  =
                                                   2     > gi+2 − gi+1
                                                         >             ,   αi+ 1 < 0.
                                                             gi+1 − gi
                                                         :
                                                                              2




                                                                                        Benasque09 – p.23/47
A Well Balanced second order Hybrid scheme
                                                        △t                  △t n
Second order method                  Uin+1 = Uin −      △x
                                                             n
                                                           (Gi+ 1  n
                                                               − Gi− 1 ) − △x Si
                                                                2    2
                                                                                    ˛n
                                 “                              ”                ∂f ˛
    n       HI
   Gi+ 1 = Gi+ 1 =           1
                             2
                                    n      n
                                                i+ 2
                                                      n      n
                                   gi+1 + gi − αn 1 (gi+1 − gi ) ,      αn 1 = λ ˛
                                                                          i+ 2
       2             2                                                           ∂u ˛i+ 1
                                                                                         2



                                  n       LO             HI      LO
                                 Gi+ 1 = Gi+ 1 + φi+ 1 (Gi+ 1 − Gi+ 1 )
                                     2         2        2      2     2

                    “                                     ”
     LO         1     n      n          n      n      n
    Gi+ 1   =   2
                     gi+1 + gi − sign(αi+ 1 )(gi+1 − gi )
        2                                  2
                    “                               ”
     HI         1     n      n     n     n      n
    Gi+ 1 =     2
                     gi+1 + gi − αi+ 1 (gi+1 − gi )
        2                                  2

    φi+ 1 = φ(ri+ 1 )
        2                2

                                                         > gi − gi−1 ,
                                                         8
                                                         >
                                                         < g               αi+ 1 > 0;
                                                             i+1 − gi         2
    Upwind smoothness indicator ri+ 1                  =
                                                   2     > gi+2 − gi+1
                                                         >             ,   αi+ 1 < 0.
                                                             gi+1 − gi
                                                         :
                                                                              2



                         gi+1 = gi ,     ∀i → U n+1 = U n Well Balanced

                                                                                        Benasque09 – p.23/47
Greenberg&Leroux tests [SINUM,96]
                                                                   (
                               2
                              u                                              1.3,       x<0
                      ut +      = −ax (x)u,             u+a=
                              2                                              1          x>0

                                   the WB second order scheme
 1.4                               1.4                                 1.5




 1.2                               1.2




  1                                 1
                                                                        1


 0.8                               0.8




 0.6                               0.6


                                                                       0.5
 0.4                               0.4




 0.2                               0.2




  0                                 0                                   0
 −0.5   0   0.5   1     1.5    2   −0.5   0   0.5   1    1.5   2       −0.5         0   0.5   1   1.5   2


 1.4                               1.4                                 1.4




 1.2                               1.2                                 1.2




  1                                 1                                   1




 0.8                               0.8                                 0.8




 0.6                               0.6                                 0.6




 0.4                               0.4                                 0.4




 0.2                               0.2                                 0.2




  0                                 0                                   0
 −0.5   0   0.5   1     1.5    2   −0.5   0   0.5   1    1.5   2       −0.5         0   0.5   1   1.5   2




            Curbing the oscillations: the WB second order hybrid scheme

                                                                                                            Benasque09 – p.24/47
Embid’s Problem
Scalar model for gas flow through a duct of variable cross section.
                               (
                                     wt + w2 /2 x = (6x − 3)w,
                                           `     ´
                                                                                                    0<x<1
                                     w(0, t) = 1; w(1, t) = −.1

                                               (
                                                     1 + 3x2 − 3x,                       x < xj
Steady solution: w(x) =
                                                     −0.1 + 3x2 − 3x,                    x > xj

        1                                                                     1


       0.8                                                                   0.8


       0.6                                                                   0.6


       0.4                                                                   0.4


       0.2                                                                   0.2


        0                                                                     0


      −0.2                                                                  −0.2


      −0.4                                                                  −0.4


      −0.6                                                                  −0.6


      −0.8                                                                  −0.8


       −1                                                                    −1
             0   0.1   0.2   0.3   0.4   0.5   0.6    0.7   0.8   0.9   1          0   0.1   0.2   0.3   0.4   0.5   0.6   0.7   0.8   0.9   1



                             WB2-scheme                                                 Curbing the oscillations
                                                                                                                                                 Benasque09 – p.25/47
FSB for Shallow Water Flows
Incorporate the source term in the flux divergence for shallow water flows
                                   0          1 0             1
                                         0               0
Split source term, S = S1 + S2 = @ −ghzx A + @           0
                                   B          C B             C
                                                              A
                                         0             −ghzy
                                      Z       x
                     B(x, y, t) = −               S1 (s, y, t)ds;   ∂x B = −S1
                                          0
 Define functions:                     Z       y
                     C(x, y, t) = −               S2 (x, s, t)ds;   ∂y C = −S2
                                          0

Formally Rewrite the original system in ’conservation form’ as follows,

                         Ut + (F + B)x + (E + C)y = 0

                              Ut + (G)x + (H)y = 0


→ combined fluxes: physical fluxes + primitive of source term


                                                                                 Benasque09 – p.26/47
The WBH2 scheme for 1D Shallow water flows

                                        !                         !                    !
                                    h                    q                      0
  Ut + (F + B)x = 0 , U =                   ,F =   q2
                                                                      B=   Rx
                                    q              h
                                                        + 2 gh2
                                                          1
                                                                                ghzx

                                     ∂F
Characteristic-based extension:         = J = RΛL, Spectral decomp. of J
                                     ∂U

                                     △t ` n           ´ △t n
                 Uin+1   =   Uin   −     G       n
                                              − Gi−1/2 −    Si
                                     △x i+1/2            △x




                                                                                   Benasque09 – p.27/47
The WBH2 scheme for 1D Shallow water flows

                                        !                              !                    !
                                    h                         q                      0
  Ut + (F + B)x = 0 , U =                   ,F =        q2
                                                                           B=   Rx
                                    q                   h
                                                             + 2 gh2
                                                               1
                                                                                     ghzx

                                     ∂F
Characteristic-based extension:         = J = RΛL, Spectral decomp. of J
                                     ∂U

                                     △t ` n           ´ △t n
                 Uin+1   =   Uin   −     G       n
                                              − Gi−1/2 −    Si
                                     △x i+1/2            △x
                                            2
                                            X
                               n                   n,p    n,p
                              Gi+1/2 =            Gi+1/2 Ri+1/2 ,
                                            p=1

                 ∗
     Ji+1/2 = J(Ui+1/2 ), U ∗ Roe’s average state at the interface.
     Gn,p = Ln,p Gn , αn 1 = p-th eigenvalue of Ji+1/2 .
      j      i+1/2 j   i+     2
      n,p
     Gi+1/2 : WBH2 p-th characteristic flux, constructed ’as in the scalar case’



                                                                                        Benasque09 – p.27/47
The WBH2 scheme for 1D Shallow water flows

                                          !                         !                         !
                                      h                    q                           0
  Ut + (F + B)x = 0 , U =                     ,F =   q2
                                                                        B=        Rx
                                      q              h
                                                          + 2 gh2
                                                            1
                                                                                       ghzx

                                       ∂F
Characteristic-based extension:           = J = RΛL, Spectral decomp. of J
                                       ∂U

                                       △t ` n           ´ △t n
                   Uin+1   =   Uin   −     G       n
                                                − Gi−1/2 −    Si
                                       △x i+1/2            △x
                  n    1“ n     n        n    n      n       n
                                                                ”
                 Si =     β 1 (Gi+1 − Gi ) + βi− 1 (Gi − Gi−1 )
                       2 i+ 2                    2

                                              !
                   ∂S ˛n           0       0            ˛       zi+1 − zi
        βi+1/2   =    ˛ 1 =         ˛n               zx ˛i+ 1 ≈
                   ∂U i+ 2     −gzx ˛i+ 1 0                 2      ∆x
                                                 2

                                                                                   !
                                                                             0
          Gi+1 − Gi = (Fi+1 − Fin + Bi,i+1 ),
                        n            n                      n
                                                           Bi,i+1 =
                                                                         bi,i+1


                                                                                          Benasque09 – p.27/47
The C-property
         Numerical integration crucial for Well Balancing
                                                 !
                                                                        xi+1
                                          0
                                                                    Z
                         n
                        Bi,i+1   =                   ,   bi,i+1 =              ghzx dx
                                        bi,i+1                      xi


                                                      g
         Trapezoidal rule (as in GL-problems) ˆi,i+1 = (zi+1 − zi )(hi + hi+1 ),
                                              b
                                                      2

                    bi,i+1 = ˆi,i+1 + O(△x3 )
                             b                             →    Approx. C- prop.
for quiescent steady flows: q = 0, h + z =constant,
Z    xi+1               Z    xi+1
                                                  g            g
            ghzx dx =               gh(−h)x dx = − (h2 − h2 ) = (zi+1 − zi )(hi + hi+1 )
                                                     i    i+1
    xi                      xi                    2            2

                 bi,i+1 = ˆi,i+1
                          b             ⇒        the WBH2 scheme is exactly WB




                                                                                         Benasque09 – p.28/47
h+z =c             q = 0 : Exact WB

                                       2.5
                                                                                   multiresolution solution
                                                                                   topography

                                        2

           (−.4∗(x−10)2 )
z(x) = .2e                             1.5



Initial data: h + z = 2, q = 0          1



L1 -error after 50 s ≈ 3 · 10−14       0.5




                                        0
                                             0   2   4   6         8   10   12      14     16      18         20




                                       14


                                       12


Non-smooth bottom topography (tab-     10


ulated)                                 8
                                                                                   multiresolution solution
                                                                                   topography

Initial data: h + z = 12                6


                                        4
                                 −15
L1 -error after 200 s ≈ 1 · 10          2


                                        0
                                             0               500                 1000                     1500




                                                                                                                   Benasque09 – p.29/47
Transcritical flow over a hump
                                   (
                                            .2 − .05(x − 10)2 ,          8 < x < 12
                          z(x) =
                                            0,                              else


             upstream q = 1.53                                            upstream q = .18

   1.4                                                        0.45
                                   numerical solution                                      numerical solution
                                   analytical solution         0.4                         analytical solution
   1.2                             topography                                              topography

                                                              0.35
    1
                                                               0.3

   0.8                                                        0.25


   0.6                                                         0.2


                                                              0.15
   0.4
                                                               0.1

   0.2
                                                              0.05


    0                                                           0
         0   5       10      15        20                25          0    5     10    15      20                 25




                 Without shock                                                With shock




                                                                                                                      Benasque09 – p.30/47
Dam Break over a discontinuous topography

Surface level h + z. Left: Numerical h + z using 400 grid cells, and bottom topography.
Right: Numerical h + z using 400 and 4000 grid cells. Top: t = 15s, Bottom t = 60s

                                           22                                                                  21
                                                                                                                                         N=4000
                                           20                                                                                            N=400
                                                                                                               20
                                           18
              surface level h+z,bottom z




                                           16                                                                  19
                                           14




                                                                                           surface level h+z
                                                                                                               18
                                           12

                                           10
                                                                                                               17
                                            8

                                            6                                                                  16

                                            4
                                                                                                               15
                                            2

                                            0                                                                  14
                                                0   500       1000   1500                                           0   500       1000        1500
                                                          x                                                                   x




                                           22                                                                  20
                                                                                                                                         N=4000
                                           20                                                        19.5                                N=400

                                           18                                                                  19
              surface level h+z,bottom z




                                           16
                                                                                                     18.5
                                           14
                                                                            surface level h+z


                                                                                                               18
                                           12
                                                                                                     17.5
                                           10
                                                                                                               17
                                            8
                                                                                                     16.5
                                            6

                                            4                                                                  16

                                            2                                                        15.5

                                            0                                                                  15
                                                0   500       1000   1500                                           0   500       1000        1500
                                                          x                                                                   x




                                                                                                                                                     Benasque09 – p.31/47
Le Veque’s quasi-stationary flow

ǫ- perturbation of quiescent flow, ǫ = 10−3


                     1.0005
                                        1st order
                                        minmod
                                        superbee
                     1.0004



                     1.0003



                     1.0002



                     1.0001



                         1



                     0.9999
                              0   0.1      0.2      0.3   0.4   0.5   0.6   0.7   0.8   0.9   1




              first order, WBH2, minmod-limiter, WBH2, Roe-limiter



                                                                                                  Benasque09 – p.32/47
2D-extension

                          Ut + (F + B)x + (E + C)y = 0

                               Ut + (G)x + (H)y = 0


→ combined fluxes: physical fluxes + primitive of source term

                          Gi+ 1 ,j − Gi− 1 ,j       Hi,j+ 1 − Hi,j− 1
                              2           2              2          2
               ∂t Uij +                         +                       =0
                                  ∆x                         ∆y




                                                                             Benasque09 – p.33/47
2D-extension: The C-property

                 Grid size     l1 -error
                 128 × 128   1.215 · 10−15
                 256 × 256   2.627 · 10−15
                 512 × 512   7.690 · 10−16




                                             Benasque09 – p.34/47
2-D Test: Quasi-stationary flow
                                                                         2 +(y−0.5)2 )
Same topography: z(x, y) = 0.5e−50((x−0.5)    Initial Data                                       .
                                                           
              1.01 − z(x, y), 0.1 < x < 0.2;      q1 (x, y)
   h(x, y) =                                                =0
              1 − z(x, y),    otherwise.          q2 (x, y)

      1

     0.9

     0.8                                                                                                     reference 129x129
                                                                                  1

     0.7                                                                         0.9                                                                   1.006

     0.6                                                                         0.8
                                                                                                                                                       1.005
                                                                                 0.7
     0.5                                                                                                                                               1.004
                                                                                 0.6
     0.4                                                                                                                                               1.003
                                                                                 0.5




                                                                             y
     0.3                                                                         0.4
                                                                                                                                                       1.002


     0.2                                                                         0.3                                                                   1.001

                                                                                 0.2                                                                   1
     0.1
                                                                                 0.1
                                                                                                                                                       0.999
      0
           0   0.1   0.2   0.3   0.4   0.5   0.6   0.7   0.8   0.9   1            0
                                                                                       0   0.1   0.2   0.3   0.4    0.5   0.6    0.7   0.8   0.9   1
                                                                                                                     x




                                                                                                                                                   Benasque09 – p.35/47
The Circular Dam-break Problem
[Castro et al. J. Sci. Comp. 2009]: free surface, slice of the channel at y=1. t=0,
0.05,0.15,0.25




                                                                                      Benasque09 – p.36/47
The Circular Dam-break Problem
Left: h; Right: q1 . longitudinal section at y = 1 Top: t = 0.1. Bottom: t = 0.15
800 points, 200 points

                                     0.75

                                      0.7
                                                                                                             0.2

                                     0.65
       surface level h+z, bottom z




                                                                                                             0.1
                                      0.6




                                                                                             discharge q1
                                     0.55                                                                     0


                                      0.5                                                                   −0.1

                                     0.45
                                                                                                            −0.2
                                      0.4

                                                                                                            −0.3
                                     0.35      N=200                                                                  N=200
                                               N=800                                                                  N=800
                                      0.3                                                                   −0.4
                                        0.4   0.6      0.8   1       1.2   1.4   1.6   1.8                     0.4   0.6      0.8   1       1.2   1.4   1.6   1.8
                                                                 x                                                                      x


                                     0.75

                                      0.7
                                                                                                             0.2

                                     0.65
       surface level h+z, bottom z




                                                                                                             0.1
                                      0.6                                                    discharge q1


                                     0.55                                                                     0


                                      0.5                                                                   −0.1

                                     0.45
                                                                                                            −0.2
                                      0.4

                                                                                                            −0.3
                                     0.35      N=200                                                                  N=200
                                               N=800                                                                  N=800
                                      0.3                                                                   −0.4
                                        0.4   0.6      0.8   1       1.2   1.4   1.6   1.8                     0.4   0.6      0.8   1       1.2   1.4   1.6   1.8
                                                                 x                                                                      x




                                                                                                                                                                    Benasque09 – p.37/47
[Caselles, RD, Haro Comp.& Fluids 2008]

Extending the Shu-Osher HRSC approach for FSB schemes

                       Ut + (F + B)x + (E + C)y = 0

                              Ut + (G)x + (H)y = 0

    Semi-discrete formulation on Cartesian meshes (separate spatial and
    temporal accuracy).

                            Gi+ 1 ,j − Gi− 1 ,j       Hi,j+ 1 − Hi,j− 1
                                2           2              2          2
                 ∂t Uij +                         +                       =0
                                    ∆x                         ∆y

                                           n     n+1
    Accuracy in time: Runge-Kutta schemes Uij ⇒ Uij .
    Relies on construction of Robust 1-D numerical fluxes.




                                                                               Benasque09 – p.38/47
HRSC Numerical Flux Functions
        Relies on construction of Robust 1-D numerical fluxes.
            Design first for scalar conservation laws
            High order nonlinear reconstruction of numerical fluxes for high order
            accuracy in space (ENO, WENO, PHM ...).
                                                                   ∂F
            Extend to systems via a local characteristic approach:     = J = RΛL
                                                                   ∂U
        [Shu-Osher JCP-86]
                                                  X    p     p
                                     Fi+ 1 =          Fi+ 1 Ri+ 1
                                          2               2     2
                                                  p



        [RD, Marquina JCP-96] Marquina Flux-Splitting ≡ 2-Jacobian Shu-Osher
        framework [RD, Mulet AMFM-06]

                          M
                                     X
                         Fi+ 1   =        p,+              p,−
                                         Fi+ 1 Rp (UL ) + Fi+ 1 Rp (UR )
                             2                2                     2
                                     p



 p       p,±
Fi+ 1 , Fi+ 1 characteristic-fluxes (constructed ’as in the scalar-case’)
    2       2


                                                                            Benasque09 – p.39/47
A HRSC-FSB scheme for the scalar CL
[Caselles, Haro, RD C& F 09]

Start with the scalar equation   wt + f (w)x = s(x, w),




                                                          Benasque09 – p.40/47
A HRSC-FSB scheme for the scalar CL
[Caselles, Haro, RD C& F 09]

Start with the scalar equation       wt + f (w)x = s(x, w),
                             Rx
Define            b(x, t) = − x s(y, w(y, t))dy, → bx = −s
                              ¯




                                                              Benasque09 – p.40/47
A HRSC-FSB scheme for the scalar CL
[Caselles, Haro, RD C& F 09]

Start with the scalar equation       wt + f (w)x = s(x, w),
                             Rx
Define            b(x, t) = − x s(y, w(y, t))dy, → bx = −s
                              ¯


          Rewrite      wt + gx = 0,       g(x, t) = f (w(x, t)) + b(x, t)




                                                                            Benasque09 – p.40/47
A HRSC-FSB scheme for the scalar CL
[Caselles, Haro, RD C& F 09]

Start with the scalar equation       wt + f (w)x = s(x, w),
                             Rx
Define            b(x, t) = − x s(y, w(y, t))dy, → bx = −s
                              ¯


          Rewrite      wt + gx = 0,           g(x, t) = f (w(x, t)) + b(x, t)

                                      1
Semi-discrete formulation: ∂t wi +      (gi+1/2 − gi−1/2 ) = 0
                                     ∆x




                                                                                Benasque09 – p.40/47
A HRSC-FSB scheme for the scalar CL
[Caselles, Haro, RD C& F 09]

Start with the scalar equation       wt + f (w)x = s(x, w),
                             Rx
Define            b(x, t) = − x s(y, w(y, t))dy, → bx = −s
                              ¯


           Rewrite     wt + gx = 0,           g(x, t) = f (w(x, t)) + b(x, t)

                                      1
Semi-discrete formulation: ∂t wi +      (gi+1/2 − gi−1/2 ) = 0
                                     ∆x

gi+1/2 :   Shu-Osher RF algorithm with upwind direction determined by f (w)




                                                                                Benasque09 – p.40/47
A HRSC-FSB scheme for the scalar CL
[Caselles, Haro, RD C& F 09]

Start with the scalar equation       wt + f (w)x = s(x, w),
                             Rx
Define            b(x, t) = − x s(y, w(y, t))dy, → bx = −s
                              ¯


           Rewrite        wt + gx = 0,           g(x, t) = f (w(x, t)) + b(x, t)

                                       1
Semi-discrete formulation: ∂t wi +       (gi+1/2 − gi−1/2 ) = 0
                                      ∆x

gi+1/2 :   Shu-Osher RF algorithm with upwind direction determined by f (w)
Higher order:       Use nonlinear reconstruction process on gi = fi + bi

             R xi
OBS!: bi =          s(w)dx has to be computed by numerical integration!
                                     R xl+1
     seek to involve only bl,l+1 =       xl
                                              s(w)dx



                                                                                   Benasque09 – p.40/47
[Caselles, Haro, RD C & F 09]

Shu-Osher RF algorithm:

        gi+1/2 − gi−1/2 = [ Gi+1/2 + HOTi+1/2 ] − [ Gi−1/2 + HOTi−1/2 ]




                                                                          Benasque09 – p.41/47
[Caselles, Haro, RD C & F 09]

Shu-Osher RF algorithm:

         gi+1/2 − gi−1/2 = [ Gi+1/2 + HOTi+1/2 ] − [ Gi−1/2 + HOTi−1/2 ]

Gi±1/2 first order contributions
                        8
                        > gi = fi + bi
                        >
                        <                      if f ′ > 0 in [wi , wi+1 ]
             Gi+1/2   =   gi+1 = fi+1 + bi+1   if f ′ < 0 in [wi , wi+1 ]
                        >
                        > 1 +
                        : (g + g − )           else
                           2   i   i+1


                                               (
                                                     +
                  ′                                 gi    =   (gi + αi+1/2 wi )
 αi+1/2 = max{|f (w)|, w ∈ [wi , wi+1 ]},           −
                                                   gi+1   =   (gi+1 − αi+1/2 wi+1 )

HOTi±1/2 Higher order terms, based on divided diferences of g.
Divided diferences of g only depend on bi,i+1 .




                                                                               Benasque09 – p.41/47
[Caselles, Haro, RD C & F 09]

Shu-Osher RF algorithm:

         gi+1/2 − gi−1/2 = [ Gi+1/2 + HOTi+1/2 ] − [ Gi−1/2 + HOTi−1/2 ]

Gi±1/2 first order contributions
                        8
                        > gi = fi + bi
                        >
                        <                      if f ′ > 0 in [wi , wi+1 ]
             Gi+1/2   =   gi+1 = fi+1 + bi+1   if f ′ < 0 in [wi , wi+1 ]
                        >
                        > 1 +
                        : (g + g − )           else
                           2   i   i+1



     First order scheme, applied to wt + awx = s(x, w) = aw (a > 0)
                                                      Z xi               !
             1 `                ´    1
     wt =        Gi+1/2 − Gi−1/2 =        fi − fi−1 +       s(w(z, t))dz
            ∆x                      ∆x                 xi−1
     [Roe, 86]




                                                                             Benasque09 – p.41/47
[Caselles, Haro, RD C & F 09]

Shu-Osher RF algorithm:

         gi+1/2 − gi−1/2 = [ Gi+1/2 + HOTi+1/2 ] − [ Gi−1/2 + HOTi−1/2 ]

Gi±1/2 first order contributions
                         8
                         > gi = fi + bi
                         >
                         <                                      if f ′ > 0 in [wi , wi+1 ]
             Gi+1/2    =   gi+1 = fi+1 + bi+1                   if f ′ < 0 in [wi , wi+1 ]
                         >
                         > 1 +
                         : (g + g − )                           else
                              2       i     i+1



Notice that Gi+1/2 − Gi−1/2 can be written in terms of
                  Z   xl+1        Z   xl        Z    xl+1
    bl+1 − bl =              s−            s=               s(x, w(x))dx = bl,l+1 ,      l = i − 1, i
                                                    xl




                                                                                                   Benasque09 – p.41/47
                                                 +        −
[Caselles, Haro, RD C & F 09] Gi+1/2 − Gi−1/2 = Gi+1/2 − Gi−1/2

                                                            +        −
      Gi+1/2 − Gi−1/2 = (Gi+1/2 − bi ) − (Gi−1/2 − bi ) =: Gi+1/2 − Gi−1/2




                                                                             Benasque09 – p.42/47
                                                 +        −
[Caselles, Haro, RD C & F 09] Gi+1/2 − Gi−1/2 = Gi+1/2 − Gi−1/2

                                                            +        −
      Gi+1/2 − Gi−1/2 = (Gi+1/2 − bi ) − (Gi−1/2 − bi ) =: Gi+1/2 − Gi−1/2

                    8
                    > fi
                    >
                    <                            if f ′ > 0 in [wi , wi+1 ]
          +
         Gi+1/2   =   fi+1 + bi,i+1              if f ′ < 0 in [wi , wi+1 ]
                    >
                    > 1 +
                    : (f + f − ) + 1 b           else
                      2  i      i+1 2 i,i+1
                    8
                    > fi−1 − bi−1,i
                    <                            if f ′ > 0 in [wi−1 , wi ]
          −
         Gi−1/2   =   fi                         if f ′ < 0 in [wi−1 , wi ]
                    > 1 +
                    : (f + f − ) − 1 b           else
                      2  i     i−1  2 i−1,i

                         (
                             fi+   =   fi + αi+1/2 wi ,
                             fi−   =   fi − αi−1/2 wi




                                                                              Benasque09 – p.42/47
                                                 +        −
[Caselles, Haro, RD C & F 09] Gi+1/2 − Gi−1/2 = Gi+1/2 − Gi−1/2

                                                            +        −
      Gi+1/2 − Gi−1/2 = (Gi+1/2 − bi ) − (Gi−1/2 − bi ) =: Gi+1/2 − Gi−1/2

                      8
                      > fi
                      >
                      <                                  if f ′ > 0 in [wi , wi+1 ]
          +
         Gi+1/2     =   fi+1 + bi,i+1                    if f ′ < 0 in [wi , wi+1 ]
                      >
                      > 1 +
                      : (f + f − ) + 1 b                 else
                        2  i      i+1 2 i,i+1
                      8
                      > fi−1 − bi−1,i
                      <                                  if f ′ > 0 in [wi−1 , wi ]
          −
         Gi−1/2     =   fi                               if f ′ < 0 in [wi−1 , wi ]
                      > 1 +
                      : (f + f − ) − 1 b
                        2  i     i−1  2 i−1,i
                                                         else


     +            −
                                       R xi+1
    Gi+1/2   −   Gi+1/2   = bi,i+1 =    xi
                                                −s(w(x), x)dx

     ±
    Gi+1/2 split the source term contribution at the i + 1/2 interface according to
    the wind, as specified by f ′ (u).



                                                                                      Benasque09 – p.42/47
[Caselles, Haro, RD C & F 09]


                                      +                      −
        gi+1/2 − gi−1/2        =   [ Gi+1/2 + HOTi+1/2 ] − [Gi−1/2 + HOTi−1/2 ]
                               =   G+
                                    i+1/2 − G−
                                             i−1/2

                        8
                        > fi
                        >
                        <                             if f ′ > 0 in [wi , wi+1 ]
            +
           Gi+1/2     =   fi+1 + bi,i+1               if f ′ < 0 in [wi , wi+1 ]
                        >
                        > 1 +
                        : (f + f − ) + 1 b            else
                          2  i      i+1 2 i,i+1
                        8
                        > fi−1 − bi−1,i
                        <                             if f ′ > 0 in [wi , wi+1 ]
            −
           Gi−1/2     =   fi                          if f ′ < 0 in [wl , wr ]
                        > 1 +
                        : (f + f − ) − 1 b            else
                          2  i     i−1  2 i−1,i



HOTi−1/2 Higher order terms, based on divided diferences of g.
                R xl+1
only bl,l+1 =    xl
                         s(w)dx are involved !!



                                                                                   Benasque09 – p.43/47
Extension to Systems
   For f (w) nonlinear, G±
                         i+1/2 lead to an automatic upwind splitting of the
   source term contribution at the i + 1/2 cell.




                                                                              Benasque09 – p.44/47
Extension to Systems
     For f (w) nonlinear, G±
                           i+1/2 lead to an automatic upwind splitting of the
     source term contribution at the i + 1/2 cell.

Extension to 1D-systems:
                                 1 “ +        −
                                                   ”
Ut + G(U )x = 0    ⇒    ∂t Ui +     Gi+1/2 − Gi−1/2 = 0
                                ∆x

Construction of G±
                 i+1/2 : Use spectral decomposition of
                                                         ∂F
                                                         ∂U




                                                                                Benasque09 – p.44/47
Extension to Systems
     For f (w) nonlinear, G±
                           i+1/2 lead to an automatic upwind splitting of the
     source term contribution at the i + 1/2 cell.

Extension to 1D-systems:
                                 1 “ +        −
                                                   ”
Ut + G(U )x = 0    ⇒    ∂t Ui +     Gi+1/2 − Gi−1/2 = 0
                                ∆x

Construction of G±
                 i+1/2 : Use spectral decomposition of
                                                         ∂F
                                                         ∂U


     Follow the basic design structure of Marquina Flux-Splitting technique:

                           G±             (Gp,± )L Rp (U L ) + (Gp,± )R Rp (U R )
                                      P
     The 2J scheme:         i+1/2 =     p   i+1/2                i+1/2




                                                                                Benasque09 – p.44/47
Extension to Systems
     For f (w) nonlinear, G±
                           i+1/2 lead to an automatic upwind splitting of the
     source term contribution at the i + 1/2 cell.

Extension to 1D-systems:
                                  1 “ +        −
                                                    ”
Ut + G(U )x = 0    ⇒     ∂t Ui +     Gi+1/2 − Gi−1/2 = 0
                                 ∆x

Construction of G±
                 i+1/2 : Use spectral decomposition of
                                                         ∂F
                                                         ∂U


     Follow the basic design structure of Marquina Flux-Splitting technique:

                           G±             (Gp,± )L Rp (U L ) + (Gp,± )R Rp (U R )
                                      P
     The 2J scheme:         i+1/2 =     p   i+1/2                i+1/2


     If U L = U R = U ∗ (Plain Extension of Shu-Osher to HCL+source term)
                                ±       P p,±
     The 1J Scheme:           Gi+1/2 = p Gi+1/2 Rp (U ∗ )
     U ∗ ≡ Ui+1/2 interface state.

Gp,± characteristic fluxes, computed ’as in the scalar case’ for the pth field.



                                                                                Benasque09 – p.44/47
[Caselles, Haro, RD C & F 09]

Numerical integration crucial for Well Balancing

     1J-scheme exactly well balanced for quiescent flows.
     2J-scheme only approximately well balanced, order ≥ 2.


                   1.0035

                                                                                         1.0006
                    1.003
                                                                                                                       fine mesh
                                                                                                                       1st order
                                                                                         1.0005
                   1.0025                                                                                              2nd order
                                                           fine mesh
                                                           1st order                                                   3rd order
                    1.002                                                                1.0004
                                                           2nd order
       Height[m]




                                                           3rd order
                   1.0015                                                                1.0003




                                                                             Height[m]
                    1.001                                                                1.0002


                   1.0005                                                                1.0001


                       1                                                                     1


                   0.9995                                                                0.9999
                            50   100      150        200      250      300                                  0.5                    1
                                       Position[m]                                                       Position[m]




                             The 2J-scheme                                                        The 1J-2J scheme




                                                                                                                                       Benasque09 – p.45/47
[Caselles, Haro, RD C & F 09]

                        The 1J-2J scheme can manage wet/dry fronts.


              0.7                                                                  12
                                                    t=0s
                                                    t=10s
                                                    t=20s
                                                    t=100s
              0.6                                   t=1000s
                                                    topography                     10


              0.5                                                                            t=0s
                                                                                    8        t=0.05s
                                                                                             t=0.25s




                                                                      Height [m]
              0.4                                                                            t=0.45s
 Height [m]




                                                                                             t=0.65s
                                                                                    6
                                                                                             topography

              0.3

                                                                                    4

              0.2


                                                                                    2
              0.1


                                                                                    0
               0                                                                        1      12.5           25
                    1                   12.5                     25                         Position[m]
                                      Position[m]



                         Drain on a non-flat bottom                    Dry bed generation + non-flat bottom


                        Slight loss of conservation in the 2J extension to systems.


                                                                                                          Benasque09 – p.46/47
Perspectives


   Treatment of Dry areas in the WBH2 scheme.
   Combine WBH2 with Adaptive Mesh Refinement techniques for complex
   situations.




                                                                 Benasque09 – p.47/47
Perspectives


     Treatment of Dry areas in the WBH2 scheme.
     Combine WBH2 with Adaptive Mesh Refinement techniques for complex
     situations.


Quote from J. Quirk (ICASE Report 92-7)

The concept of using several computer codes to cross-check numerical results
becomes ill-founded if all codes follow the same methodology.




                                                                         Benasque09 – p.47/47
Perspectives


     Treatment of Dry areas in the WBH2 scheme.
     Combine WBH2 with Adaptive Mesh Refinement techniques for complex
     situations.


Quote from J. Quirk (ICASE Report 92-7)

The concept of using several computer codes to cross-check numerical results
becomes ill-founded if all codes follow the same methodology.




                Thanks for your attention!


                                                                         Benasque09 – p.47/47

				
DOCUMENT INFO
Shared By:
Categories:
Tags:
Stats:
views:22
posted:11/1/2011
language:English
pages:92