Conservative Methods
Godunov's First-Order Upwind Method
$$ u^{n+1}_i = u^n_i + \frac{\Delta t}{\Delta x}\left[ f_{i-1/2} - f_{i+1/2} \right] $$
- Godunov's first-order upwind method는 위의 식의 intercell numerical flux $f_{i+1/2}$를 구할 때, local Riemann problem의 solution을 활용함. 이는 현재시간 $n$에서 값들이 piece-wise constant distribution이라는 것을 바탕으로 함.
- Intecell boundary $x_{i+1/2}$에서 discontinuity에 의해 분리되어 있는 constant states $(u^n_i, u^n_{i+1})$를 적용하여, 아래의 Riemann problem을 정의할 수 있음. \[ \begin{aligned} \text{PDE:} & \quad u_t + f(u)_x = 0, \\[5pt] \text{IC:} & \quad u(x, 0) = u_0(x) = \begin{cases} u^n_i, & x < 0, \\ u^n_{i+1}, & x > 0. \end{cases} \end{aligned} \]
First version of Godunov's method
- Cell value $u^n_i$을 새로운 value $u^{n+1}_i$로 update하기 위해서 두개의 Riemann problems, $RP(u^n_{i-1}, u^n_i)$와 $RP(u^n_i, u^n_{i+1})$,를 해결해야 함. 두개의 local problem의 solution을 integral average하여 $u^{n+1}_i$를 얻을 수 있음.
- $ RP(u^n_{i-1}, u^n_i) $ 과 $ RP(u^n_i, u^n_{i+1}) $ 의 exact solution은 아래와 같으며, 각 solution은 $(0,0)$을 원점으로 하는 local frame of reference를 가짐. \[ u_{i-\frac{1}{2}}\left(x/t \right) = \begin{cases} u^n_{i-1}, & x/t < a, \\ u^n_i, & x/t > a. \end{cases} \] \[ u_{i+\frac{1}{2}}\left(x/t \right) = \begin{cases} u^n_{i}, & x/t < a, \\ u^n_{i+1}, & x/t > a. \end{cases} \]
- Godunov scheme에 따라 *updated solution은 다음과 같음. $$ u^{n+1}_i = \frac{1}{\Delta x}\left[ \int^{\frac{1}{2}\Delta x}_0 u_{i-1/2}(x/\Delta t)dx + \int^{0}_{-\frac{1}{2}\Delta x}u_{i+1/2}(x/\Delta t) dx \right] $$
* $u_{i - 1/2}$의 right half를 쓰는 이유는 intercell $i-1/2$을 $(0,0)$으로 생각할 때 오른쪽 영역이며,
$u_{i + 1/2}$의 left half를 쓰는 이유는 intecell $i+1/2$를 기준으로 왼쪽 영역이기 때문.
- Time step에 대한 제약조건 $c=a\Delta t / \Delta x \le 1/2$를 적용하여, 다음과 같이 표현할 수 있음. 정리한 식은 CIR first order upwind method와 동일함. $$ u^{n+1}_i = u^n_i + c (u^n_{i-1} + u^n_i) $$
Second version of Godunov's method
- 또 다른 Godunov's method를 통해서, 앞서 언급되었던 over-restrictive CFL-like condition을 피할 수 있음.
- Integral average는 $RP(u^n_{i-1}, u^n_i)$와 $RP(u^n_i, u^n_{i+1})$의 combined solution인 $\tilde{u}$를 적용하여 아래와 같이 표현할 수 있음. $$ u^{n+1}_i = \frac{1}{\Delta x}\int^{i+\frac{1}{2}}_{i-\frac{1}{2}} \tilde{u}(x, \Delta t) dx $$
- Control volume $\left[ x_{i-1/2}, x_{i+1/2} \right] \times [0, \Delta t]$에서 다음과 같이 정리할 수 있음. \begin{align*}
\int^{x_{i+1/2}}_{x_{i-1/2}} \tilde{u}(x, \Delta t) \, dx &= \int^{x_{i+1/2}}_{x_{i-1/2}} \tilde{u}(x, 0) \, dx \\
&\quad + \int^{\Delta t}_{0} f(\tilde{u}(x_{i-1/2}, t)) \, dt
- \int^{\Delta t}_{0} f(\tilde{u}(x_{i+1/2}, t)) \, dt
\end{align*} - 아래의 cell average 정의와 time integral average로 정의된 *intecell flux를 활용하여 conservative formula로 만들 수 있음. $$ u^n_i = \frac{1}{\Delta x} \int^{x_{i+1/2}}_{x_{i-1/2}} u (x, t^n)dx $$ $$ f_{i\pm 1/2} = \frac{1}{\Delta t} \int^{\Delta t}_0 f(\tilde{u}(x_{i \pm 1/2}, t)) dt $$ $$ u^{n+1}_i = u^n_i + \frac{\Delta t}{\Delta x}\left[ f_{i-1/2} - f_{i+1/2} \right] $$
- 각 cell interface에서의 $f(\tilde{u}(x,t))$는 exact solution $\tilde{u}(x,t)$에 의존하며, 이는 다음과 같이 주어짐. $$ \tilde{u}(x_{x \pm 1/2}, t) = u_{i \pm 1/2}(0) $$ $$ f_{i \pm 1/2} = f(u_{i \pm 1/2} (0) ) $$
- 따라서, Godunov intercell numerical flux는 다음과 같이 표현할 수 있음. $u_{i+1/2}(0)$은 $x/t = 0$에서의 Riemann problem의 exact solution $u_{i+1/2}(x/t)$를 나타냄. 즉, intercell boundary를 따라 얻어지는 solution임. $$ f^\text{god}_{i+1/2} = f(u_{i+1/2}(0)) $$
* Intecell boundary에서 physical flux $f(u)$의 time integral average
Godunov's Method for Burgers's Equation
- Non-linear PDE에 대한 Godunov's method 적용을 이해하기 위해 inviscid Burgers equation을 다룸. 두번 째 방법의 Godunov's method를 적용하기 위해 Riemann problem의 solution $u_{i+1/2}(x/t)$이 필요함. Shock wave $(u^n_i > u^n_{i+1})$와 rarefaction wave $(u^n_i \le u^n_{i+1})$를 포함하는 solution은 다음과 같음. \[ u_{i+1/2}(x/t) = \begin{cases} u^n_{i}, & S \ge x/t, \\ u^n_{i+1}, & S \le x/t. \end{cases} \] $$ S = 1/2(u^n_i + u^n_{i+1})$$ \[ u_{i+1/2}(x/t) = \begin{cases} u^n_{i}, & x/t \le u^n_i, \\ x/t, & u^n_i < x/t < u^n_{i+1}, \\ u^n_{i+1}, & x/t \ge u^n_{i+1}. \end{cases} \]
- Godunov's flux를 구하기 위해서 $u_{i+1/2}(0)$이 필요함. 이는 local Riemann problem에서 $x/t=0$에 해당하는 intercell boundary $x_{i+1/2}$에서의 solution $u_{i+1/2}(x/t)$을 의미함.
- 또한, solution에서의 모든 가능한 wave patterns을 확인해야 함. Burgers's equation에 대해서 5개의 가능성이 존재함.
- Shock wave:
- Shock to the left
- Shock to the right
- Rarefaction wave
- Supersonic to the left
- Supersonic to the right
- Transonic rarefaction (sonic rarefaction): wave가 교차할 때, characteristic speed $u$의 부호가 변하기 때문에 $u=0$이 되는 sonic point가 존재함.
- Shock wave:
- 따라서, complete solution은 다음과 같이 정리할 수 있음. \[ u_{i+1/2}(0) = \begin{cases} u^n_{i}, & S \ge 0, \\ u^n_{i+1}, & S \le 0. \end{cases} \] $$ S = 1/2(u^n_i + u^n_{i+1})$$ \[ u_{i+1/2}(0) = \begin{cases} u^n_{i}, & 0 \le u^n_i, \\ 0, & u^n_i < 0 < u^n_{i+1}, \\ u^n_{i+1}, & 0 \ge u^n_{i+1}. \end{cases} \]
- Godunov's scheme을 Burgers's equation에 적용할 때, boundary condition를 고려해야 하며 적절한 time step $\Delta t$을 선택해야 함.
Boundary Condition
- $i=2, \dots, M-1$의 모든 cell $i$에서는 두개의 intercell fluxes가 정의되지만, left와 right boundary에 인접한 cell $1$과 $M$에서는 하나의 intercell flux를 가짐. 이를 위해서 다음의 방법들을 고려할 수 있음.
- Left boundary $x=0$에 할당되는 boundary function $u_l(t)$를 가정하여, boundary에서의 intercell flux $f_{1/2}=f(u_l(t))$를 정의함.
- Boundary $x=0$의 left에 cell average $u_0^n$을 가지는 fictitious cell $0$를 명시하여 Riemann problem $RP(u_0^n, u_1^n)$을 풀고, 이를 통해서 missing intercell flux $f_{1/2}$를 찾음.
- e.g., 아래와 같은 boundary condition을 적용하면, boundary에서 어떠한 disturbance도 생기지 않으며, boundary가 없는 것처럼, wave가 그대로 통과함. (transparent or transmissive boundary condition) $$ u^n_0 = u^n_1, \quad u^n_{M+1} = u^n_M $$
※ Fictitious states는 풀고자 하는 문제의 물리적 상황에만 의존
Time Step
- Time step $\Delta t$를 선택하는 것은 특정 scheme의 stability condition과 매우 밀접한 관련이 있음. Godunov's method의 경우, $\Delta t$의 선택은 Courant number $c$의 제약조건에 의존함.
- Non-linear problem은 각 time level에서 여러개의 wave speed가 존재하기 때문에 여러 개의 Cournat number가 존재함. 가장 빠른 wave은 time step $\Delta t$동안 최대 one cell length $\Delta x$만큼 통과함.
- Time level $n$에서 domain을 걸쳐서 지나가는 가장 빠른 wave speed를 $S^n_\text{max}$로 표시하면, 최대 Courant number $C_\text{cfl}$은 다음과 같이 정의할 수 있음. $$ 0 < C_\text{cfl} = S^n_\text{max}\frac{\Delta t}{\Delta x} \le 1, \quad \Delta t = C_\text{cfl}\frac{\Delta x}{S^n_\text{max}} $$
- Maximum speed $S^n_\text{max}$를 추정하기 위해 다음의 방법들을 고려할 수 있음.
- Fictitious cell을 포함하여 전체 $M+2$개의 characteristic speed $u_i^n$ 중에서, 절댓값으로 최대 speed를 $S^n_\text{max}$로 정함.
- 각 cell interface에서 Riemann problem으로부터 speed $S^n_{n+1/2}$를 얻고 그 중에서 최대 wave speed를 정함. Shock의 경우 shock speed를, rarefaction의 경우 expansion의 head와 tail의 speed 중에서 큰 값을 사용함. \[ S^n_{i+1/2} = \begin{cases} | \frac{1}{2}(u^n_i + u^n_{i+1}) |, & \text{shock}, \\ \text{max}(|u^n_i|, |u^n_{i+1}|), & \text{rarefaction}. \end{cases} \] \[
S^n_\text{max} = \max_i\left( S^n_{i+1/2} \right)
\]
Conservative Form of Difference Schemes
- Conventional finite difference scheme은 conservation form으로 표현될 수 있으며, 아래의 scheme을 통해서 intercell flux를 얻을 수 있음.
Lax and Friedrichs (LF) Scheme
$$ f_{i+1/2} = \frac{(1+c)}{2c} u^n_i + \frac{c-1}{2c} u^n_{i+1} $$
- LF scheme은 cell $i$안에서의 integral average로도 표현할 수 있음. $$ u^{n+1}_{i} = \frac{1}{\Delta x} \int^{x_{i+1/2}}_{ x_{i-1/2} } \tilde{u}\left( x, \frac{1}{2}\Delta t \right) dx $$ \[ \tilde{u}(x/t) = \begin{cases} u^n_{i-1}, & x/t < a, \\ u^n_{i+1}, & x/t > a. \end{cases} \]
- LF solution은 Riemann problem solution의 weighted average이며, upwind term이 항상 larger weight이기 때문에 upwind biased임.
- Riemann problem $RP(\mathbf{U}^{n}_{i-1}, \mathbf{U}^{n}_{i+1} )$의 solution $\tilde{\mathbf{U}}(x,t)$를 활용하여, conservation laws의 non-linear system에 LF를 적용하면 위의 식은 아래와 같이 됨. $$ \mathbf{U}^{n+1}_i = \frac{1}{\Delta x}\int^{x_{i+1/2}}_{x_{i-1/2}} \tilde{\mathbf{U}}\left( x, \frac{1}{2}\Delta x \right) dx $$ \[
\begin{aligned}
\int^{\frac{1}{2}\Delta x}_{-\frac{1}{2}\Delta x} \tilde{\mathbf{U}}\left( x, \frac{1}{2}\Delta x \right) dx
&= \int^{\frac{1}{2}\Delta x}_{-\frac{1}{2}\Delta x} \tilde{\mathbf{U}}\left( x, 0 \right) dx \\
&+ \int^{\frac{1}{2}\Delta t}_0 \mathbf{F}\left(\tilde{\mathbf{U}}\left(-\frac{1}{2}\Delta x, t\right)\right) dt \\
&- \int^{\frac{1}{2}\Delta t}_0 \mathbf{F}\left(\tilde{\mathbf{U}}\left(\frac{1}{2}\Delta x, t\right)\right) dt
\end{aligned}
\] $$ \mathbf{U}^{n+1}_i = \frac{1}{2}\left( \mathbf{U}^{n}_{i-1} + \mathbf{U}^{n}_{i+1} \right) + \frac{1\Delta t}{2\Delta x} \left( \mathbf{F}^{n}_{i-1} - \mathbf{F}^{n}_{i+1} \right) $$ - 따라서, LF scheme의 conservative version으로 정리하면 다음과 같음. $$ \mathbf{U}^{n+1}_{i} = \mathbf{U}^{n}_{i} + \frac{\Delta t}{\Delta x} \left[ \mathbf{F}_{i-1/2} - \mathbf{F}_{i+1/2} \right] $$ $$ \mathbf{F}^\text{LF}_{i+1/2} = \frac{1}{2}(\mathbf{F}^n_i + \mathbf{F}^n_{i+1}) + \frac{1\Delta x}{2\Delta t}(\mathbf{U}^n_i - \mathbf{U}^n_{i+1}) $$
Lax and Wendroff (LW) Scheme
$$ f_{i+1/2} = \frac{(1+c)}{2} au^n_i + \frac{1-c}{2} au^n_{i+1} $$
- Linear advection equation에 대해서 LW scheme은 다음과 같이 표현할 수 있음. $$ f_{i+1/2} = f(u^{n+1/2}_{i+1/2}); \quad u^{n+1/2}_{i+1/2}= \frac{(1+c)}{2} u^n_i + \frac{1-c}{2} u^n_{i+1} $$
- Burgers's equation과 같은 non-linear scalar conservation laws에 대해서는 다음과 같이 일반화할 수 있음. $$ f_{i+1/2} = f(u^{n+1/2}_{i+1/2}); \quad u^{n+1/2}_{i+1/2} = \frac{1}{\Delta x}\int^{1/2\Delta x}_{-1/2\Delta x} u_{i+1/2}\left( x, \frac{1}{2}\Delta t \right) dx $$
- Non-linear hyperbolic system에 대해서는 다음과 같이 쓸 수 있으며, 이를 Weighted Average Flux (WAF) method라고 함. $$ \mathbf{F}_{i+1/2} = \mathbf{F}\left( \mathbf{U}^{n+1/2}_{i+1/2} \right); \quad \mathbf{U}^{n+1/2}_{i+1/2}=\frac{1}{\Delta x}\int^{1/2\Delta x}_{-1/2\Delta x}\mathbf{U}_{i+1/2}\left( x, \frac{1}{2}\Delta t \right) dx$$ $$ \mathbf{F}_{i+1/2} = \mathbf{F}\left( \mathbf{U}^{n+1/2}_{i+1/2} \right); \quad \mathbf{U}^{n+1/2}_{i+1/2} = \frac{1}{2}\left( \mathbf{U}^{n}_{i} + \mathbf{U}^{n}_{i+1} \right) + \frac{1\Delta t}{2\Delta x} \left( \mathbf{F}^{n}_{i} - \mathbf{F}^{n}_{i+1} \right) $$
Upwind Schemes for Linear Systems
- Constant coefficient를 가지는 hyperbolic system은 characteristic variable로 표현되는 decoupled system으로 나타낼 수 있음. $$ \mathbf{U}_t + \mathbf{AU}_x = 0 $$ $$ \mathbf{V}_t + \mathbf{\Lambda V}_x = 0 $$ $$ \frac{\partial}{\partial t}v_j + \lambda_j \frac{\partial}{\partial x}v_j = 0 $$
- One-sided differencing scheme은 모든 eigenvalues가 같은 부호일 때만 유효함. 하지만 일반적으로 eigenvalues의 부호는 혼합되어 있기 때문에, upwind와 downwind side를 적절히 선택해야 함. 이는 eigenvalue matrix를 positivie 또는 zero eigenvalue를 가지는 matrix와 negative 또는 zero eigenvalue를 가지는 matrix로 나눠서 해결할 수 있음.
CIR Scheme
- CIR scheme을 decoupled linear hyperbolic system에 적용하면 다음과 같음. \[
\begin{aligned}
(v_j)^{n+1}_i = (v_j)^n_i &- \frac{\Delta t}{\Delta x}\lambda^+_j \left[(v_j)^n_i - (v_j)^n_{i-1}\right] \\
&\quad - \frac{\Delta t}{\Delta x}\lambda^-_j \left[(v_j)^n_{i+1} - (v_j)^n_{i}\right]
\end{aligned}
\] \[
\begin{aligned}
\mathbf{V}^{n+1}_i = \mathbf{V}^n_i &- \frac{\Delta t}{\Delta x}\mathbf{ \Lambda }^+ \left(\mathbf{V}^n_i - \mathbf{V}^n_{i-1}\right) \\
&\quad - \frac{\Delta t}{\Delta x}\mathbf{ \Lambda }^- \left(\mathbf{V}^n_{i+1} - \mathbf{V}^n_{i}\right)
\end{aligned}
\] - Original variables $\mathbf{U} = \mathbf{KV}$로 표현하면 다음과 같음. \[
\begin{aligned}
\mathbf{U}^{n+1}_i = \mathbf{U}^n_i &- \frac{\Delta t}{\Delta x}\mathbf{A}^+ \left(\mathbf{U}^n_i - \mathbf{U}^n_{i-1}\right) \\
&\quad - \frac{\Delta t}{\Delta x}\mathbf{ A }^- \left(\mathbf{U}^n_{i+1} - \mathbf{U}^n_{i}\right)
\end{aligned}
\] - Coefficient matrix $\mathbf{A}$를 positive와 negative component로 분리함으로써, characteristic speed의 부호에 따라 one-sided spatial differencing을 적용할 수 있음.
Godunov's Method
- Godunov's first-order upwind method를 통해 linear hyperbolic system을 풀기 위해서는 아래와 같이 intercell numerical flux를 계산해야함. 이는 local Riemann problem $RP(\mathbf{U}^n_i, \mathbf{U}^n_{i+1})$ 의 solution $\mathbf{U}_{i+1/2}(x/t)$이 필요함. $$ \mathbf{F}_{i+1/2} = \mathbf{F}(\mathbf{U}_{i+1/2}(0)) $$
- Solution $\mathbf{U}_{i+1/2}(x/t)$는 initial data $\mathbf{U}^{n}_i, \mathbf{U}^n_{i+1}$를 통해 다음과 같이 얻을 수 있음. $$ \mathbf{U}^n_i = \sum^m_{j=1}\alpha_j \mathbf{K}^{(j)}, \quad \mathbf{U}^n_{i+1} = \sum^m_{j=1}\beta_j\mathbf{K}^{(j)} $$ $$ \mathbf{U}_{j+1/2}(x/t) = \sum^I_{j=1}\beta_j\mathbf{K}^{(j)} + \sum^m_{j=I+1}\alpha_j\mathbf{K}^{(j)} $$
- Godunov flux는 $x/t=0$에서의 solution이 필요함. $x/t=0$의 위치는 $\lambda_I \le 0$과 $\lambda_{I+1} \ge 0$ 사이이기 때문에 다음과 같이 solution $\mathbf{U}_{i+1/2}(0)$을 얻을 수 있음. $$ \mathbf{U}_{i+1/2}(0) = \mathbf{U}^{n}_i + \sum^I_{j=1} (\beta_j - \alpha_j)\mathbf{K}^{(j)} $$
- 위의 식은 left data state $\mathbf{U}_i^n$에 negative 또는 zero speed를 가지는 wave에 걸친 모든 wave jumps를 더하면 $x/t=0$에서 Riemann problem의 solution이라는 것을 보여줌. $$ \mathbf{U}_{i+1/2}(0) = \mathbf{U}^{n}_{i+1} - \sum^m_{j=I+1} (\beta_j - \alpha_j)\mathbf{K}^{(j)} $$
- 유사하게, 위의 식은 right data state $\mathbf{U}_{i+1}^n$에 positive 또는 zero speed를 가지는 wave에 걸친 모든 wave jumps를 빼면 solution이 된다는 것을 보여줌. 위의 두 식을 합치면 다음과 같이 표현할 수 있음. $$ \mathbf{U}^{i+1/2}(0) = \frac{1}{2}\left( \mathbf{U}^n_i + \mathbf{U}^n_{i+1} \right) - \frac{1}{2} \sum^m_{j=1}\text{sign}(\lambda_j)(\beta_j - \alpha_j)\mathbf{K}^{(j)} $$
- 위의 solution을 활용하여, Godunov intercell numerical flux $\mathbf{F(U)}$를 다음과 같이 얻을 수 있음. $$ \mathbf{F}_{i+1/2} = \mathbf{F}^n_i + \sum^I_{j=1}\mathbf{A}(\beta_j - \alpha_j)\mathbf{K}^{(j)} $$ $$ \mathbf{F}_{i+1/2} = \mathbf{F}^n_i + \sum^I_{j=1}(\beta_j - \alpha_j)\lambda_j\mathbf{K}^{(j)} $$ $$ \mathbf{F}_{i+1/2} = \mathbf{F}^n_i - \sum^{m}_{j=I+1}(\beta_j - \alpha_j)\lambda_j\mathbf{K}^{(j)} $$ $$ \mathbf{F}_{i+1/2} = \frac{1}{2}(\mathbf{F}^n_i + \mathbf{F}^n_{i+1}) - \frac{1}{2}\sum^m_{j=1}(\beta_j - \alpha_j)|\lambda_j|\mathbf{K}^{(j)} $$
- 또한, Godunov flux는 또 다른 형태로 표현할 수 있음. $$ \mathbf{F}_{i+1/2} = \frac{1}{2}(\mathbf{F}^n_i + \mathbf{F}^n_{i+1}) - \frac{1}{2}|\mathbf{A}|(\mathbf{U}^n_{i+1} - \mathbf{U}^n_i) $$ \[
\mathbf{F}_{i+1/2} = \mathbf{A}^+\mathbf{U}^n_i + \mathbf{A}^-\mathbf{U}^n_{i+1}
\]