논문 링크:
https://www.sciencedirect.com/science/article/abs/pii/S002199912300462X?via%3Dihub
논문에서는 푸아송 방정식과 같은 타원 문제를 불연속 갤러킨(DG) 방법으로 해석할 때, 기울기(gradient)의 수렴성을 한 차원 개선하는 새로운 접근법을 제시합니다. 이 방법은 기존에 하이퍼볼릭 재구성(hyperbolic reformulation)을 통해 얻은 장점을 DG 방법에 접목시키고, 이를 Summation‐by‐Parts (SBP) operator를 이용하여 국소적으로 기울기를 계산하는 방식으로 구현합니다. 아래에 수식을 최대한 반영하여 상세하게 설명합니다.
1. 하이퍼볼릭 재구성과 문제의 약식 표현
기본적으로 푸아송 방정식
$$
-\Delta \phi = f
$$
를 풀기 위해, Nishikawa [17]는 이 타원 문제를 아래와 같은 선형 하이퍼볼릭 시스템의 정상 상태(steady state) 문제로 재구성합니다:
$$
\begin{aligned}
\phi_t - \nabla \cdot \mathbf{q} &= f, \\
\mathbf{q}_t - \frac{1}{\tau}\nabla \phi &= -\frac{1}{\tau}\mathbf{q},
\end{aligned}
\tag{1.1}
$$
여기서 $\tau>0$는 수렴 속도를 조절할 수 있는 완화(relaxation) 시간입니다. 정상 상태에 도달하면 시간 미분 항이 소멸하여, DG 방법의 공간 이산화에 초점을 맞출 수 있습니다.
2. 약한 형식(Weak Formulation) 및 수치 플럭스
하이퍼볼릭 시스템의 정상 상태를 대상으로 DG 방법을 적용하면, 각 요소(element) $I=(x_G,x_{G+1})$에서 시험함수 $\varphi,\psi$를 사용한 약한 형식은 다음과 같이 쓸 수 있습니다:
$$
\int_{x_G}^{x_{G+1}} \mathbf{q}\, \varphi \,dx = \int_{x_G}^{x_{G+1}} f\, \varphi \,dx + \Bigl[\,\widehat{\mathbf{q}}\, \varphi\,\Bigr]_{x_G}^{x_{G+1}},
$$
$$
\int_{x_G}^{x_{G+1}} \mathbf{q}\, \psi \,dx = -\int_{x_G}^{x_{G+1}} \phi\, \psi' \,dx + \Bigl[\,\widehat{\phi}\, \psi\,\Bigr]_{x_G}^{x_{G+1}},
\tag{2.1}
$$
여기서 $\widehat{\phi}$와 $\widehat{\mathbf{q}}$는 요소 경계에서 정의되는 수치 플럭스(numerical flux)입니다.
수치 플럭스로는 하이퍼볼릭 시스템에 자연스럽게 등장하는 upwind 플럭스를 사용하며, 이는
$$
\widehat{\phi} = \{\!\!\{\phi\}\!\!\} + \frac{\sqrt{\tau}}{2}\, [[\mathbf{q}]], \quad
\widehat{\mathbf{q}} = \{\!\!\{\mathbf{q}\}\!\!\} + \frac{1}{2\sqrt{\tau}}\, [[\phi]],
\tag{2.2}
$$
로 정의됩니다.
여기서
- $\{\!\!\{\cdot\}\!\!\}$는 요소 양쪽의 값의 산술 평균(average)을,
- $[[\cdot]]$는 요소 경계에서의 값의 차이(점프, jump)를 나타냅니다.
또한 Castillo et al. [5]에서 사용된 보다 일반적인 플럭스 형식
$$
\widehat{\phi} = \{\!\!\{\phi\}\!\!\} - \alpha_{12}\, [[\phi]] + \alpha_{22}\, [[\mathbf{q}]], \quad
\widehat{\mathbf{q}} = \{\!\!\{\mathbf{q}\}\!\!\} + \alpha_{11}\, [[\phi]] + \alpha_{12}\, [[\mathbf{q}]]
\tag{2.3}
$$
와 비교하면, (2.2)는 $\alpha_{12}=0$인 특수한 경우임을 알 수 있습니다. Cockburn, Guzmán, Wang [7]는 이 설정 하에서 DG 방법이 최적의 오차 추정(optimal error estimates)을 보인다고 분석하였습니다.
3. SBP Operators와 DG Discretization
각 요소 $\kappa$ 내에서 DG 방법은 다음과 같은 구성요소로 이루어집니다:
1. 미분 근사자 $D_\kappa$:
연속 미분 연산자 $\partial_x$를 근사합니다.
2. 대각 질량 행렬 $M_\kappa$:
$L^2$ 내적을 근사하는 역할을 하며, 대각 성질 때문에 경계 노드에서의 처리가 용이합니다.
3. 경계 보간 연산자 $C_{\kappa,L}$와 $C_{\kappa,R}$:
요소의 왼쪽/오른쪽 경계에서 다항식 근사값을 평가합니다. 예를 들어, DG spectral element method(DGSEM)의 경우, Gauss–Lobatto–Legendre (GLL) 노드를 사용하므로 경계 노드가 포함됩니다.
이들 SBP operator는 “integration by parts”를 모방하는 SBP 조건
$$
M_\kappa\, D_\kappa + \bigl(M_\kappa\, D_\kappa\bigr)^T = C_{\kappa,R}^T\, C_{\kappa,R} - C_{\kappa,L}^T\, C_{\kappa,L}
$$
를 만족합니다. 이 조건 덕분에 경계에서 발생하는 표면 항(surface terms)이 국소적으로 명시적으로 취급될 수 있습니다.
DG 이산화의 강형식(strong form)은 (2.1)을 요소 내에서 다시 기술하면
$$
\mathbf{q}_\kappa = D_\kappa\, \phi_\kappa + M_\kappa^{-1} C_{\kappa,R}\Bigl( \widehat{\phi} - C_{\kappa,R}\, \phi_\kappa \Bigr) - M_\kappa^{-1} C_{\kappa,L}\Bigl( \widehat{\phi} - C_{\kappa,L}\, \phi_\kappa \Bigr),
$$
$$
-\, D_\kappa\, \mathbf{q}_\kappa = f + M_\kappa^{-1} C_{\kappa,R}\Bigl( \widehat{\mathbf{q}} - C_{\kappa,R}\, \mathbf{q}_\kappa \Bigr) - M_\kappa^{-1} C_{\kappa,L}\Bigl( \widehat{\mathbf{q}} - C_{\kappa,L}\, \mathbf{q}_\kappa \Bigr).
\tag{2.4}
$$
이때, $\phi_\kappa$와 $\mathbf{q}_\kappa$는 요소 $\kappa$ 내에서의 다항식 근사를 나타내는 계수 벡터입니다.
4. 기울기(@)의 국소 계산 및 Theorem 2.1
논문의 핵심 정리(Theorem 2.1)는 대각 노름(diagonal-norm) SBP operator를 사용할 경우,
인접 요소의 경계에서의 해와 기울기 값만으로 내부 요소의 기울기 $\mathbf{q}$를 국소적으로 계산할 수 있음을 보입니다.
4.1. 두 요소를 통한 분석
두 요소를 좌측(;)와 우측(A)으로 나누어 고려하면, 각 요소에서의 강형식 표현은 (2.5)와 같이 나타납니다.
예를 들어, 좌측 요소에서는
$$
\mathbf{q}_{;}= D_{;} \phi_{;} + M_{;}^{-1} C_{;,R}\Bigl(\widehat{\phi} - C_{;,R}\,\phi_{;}\Bigr) + \text{추가 표면 항},
\tag{2.5}
$$
그리고 우측 요소에서는 유사한 형태의 식이 주어집니다.
경계에서의 값(예, $C_{;,R}\phi_{;}$와 $C_{A,L}\phi_{A}$)은 경계 노드에서의 다항식 근사값을 나타내며, 대각 질량 행렬의 특성 때문에 요소 내부의 다른 표면 항은 소멸합니다.
두 요소의 경계에서의 기울기 값의 차이(점프)는 다음과 같이 표현됩니다:
$$
C_{A,L}\, \mathbf{q}_{A} - C_{;,R}\, \mathbf{q}_{;} = -\frac{1}{2}\beta_{1} \Bigl( C_{A,R}\,\phi_{A} - C_{;,R}\,\phi_{;}\Bigr) + \frac{1}{2}\beta_{2}\Bigl( C_{A,L}\,\mathbf{q}_{A} - C_{;,R}\,\mathbf{q}_{;}\Bigr),
\tag{2.7}
$$
여기서 $\beta_{1}$와 $\beta_{2}$는 SBP operator의 구성(특히 그 격자(grid)와 보간 연산자의 특성)에 따라 결정되며, 균일 격자에서는 $\beta_{1}$가 0가 되는 특성이 있습니다.
이 관계를 이용하면, 각 요소 내에서 기울기 $\mathbf{q}_\kappa$를 다음과 같이 명시적으로 계산할 수 있습니다:
$$
\begin{aligned}
\mathbf{q}_\kappa =\; & D_\kappa\, \phi_\kappa + M_\kappa^{-1} C_{\kappa,R}\left(\frac{1}{2}(1-\sqrt{\tau}\,\beta_{1,R})\, [[\phi]]_{R} + \frac{\sqrt{\tau}}{2}\,\beta_{2,R}\, [[\mathbf{q}]]_{R}\right)\\[1mm]
& + M_\kappa^{-1} C_{\kappa,L}\left(\frac{1}{2}(1+\sqrt{\tau}\,\beta_{1,L})\, [[\phi]]_{L} - \frac{\sqrt{\tau}}{2}\,\beta_{2,L}\, [[\mathbf{q}]]_{L}\right).
\end{aligned}
\tag{2.10}
$$
여기서
- $[[ \phi]]_{R} = C_{A,R}\,\phi_{A} - C_{;,R}\,\phi_{;}$ (오른쪽 경계에서의 점프),
- $[[ \mathbf{q}]]_{R} = C_{A,L}\,\mathbf{q}_{A} - C_{;,R}\,\mathbf{q}_{;}$ (유사하게 정의됨)
입니다.
이 식은 SBP 구조 덕분에 요소 내부에서 기울기를 인접 요소의 경계값만을 사용하여 국소적으로 계산할 수 있음을 보여주며, 전역 시스템을 풀 필요가 없게 됩니다.
4.2. 경계 조건의 처리
경계에서 Dirichlet 조건이 약하게 부과될 경우, 에너지 분석에 기반한 수치 플럭스는 다음과 같이 사용됩니다:
$$
\widehat{\phi} = \phi_{\text{boundary}}, \quad \widehat{\mathbf{q}} = \mathbf{q}_{\text{interior}} + \frac{\sqrt{\tau}}{2}\, [[\phi]].
\tag{2.11}
$$
이로 인해, 경계 요소에서는 내부 요소와 동일한 방식으로 (2.12)와 (2.13) 식에 따라 기울기를 계산할 수 있습니다.
5. 수치 실험과 수렴률
논문 후반부에서는 1D 및 2D 문제에 대해 제시한 방법의 성능을 검증합니다.
- 1D 실험:
다항식 차수 $p$를 사용했을 때, 해 $\phi$와 기울기 $\mathbf{q}$ 모두 $p+1$ 차의 수렴률을 보였습니다. (예, $p=2$인 경우 수렴률이 약 3, $p=3$인 경우 약 4)
- 2D 실험:
2D 문제에서도 동일하게 최적 수렴률이 관찰되었으며, 경계 조건과 격자 설정에 따른 효과를 확인할 수 있었습니다.
이러한 결과는 DG 방법에 SBP operator를 결합함으로써, 기존에 기울기 계산에 있어서 전역 선형 시스템을 풀어야 하는 어려움을 해소하고, 기울기의 superconvergence 즉, 한 차수 높은 수렴률을 달성할 수 있음을 보여줍니다.
결론
요약하면, 본 논문은
1. 푸아송 방정식을 하이퍼볼릭 시스템으로 재구성하여 DG 방법에 적용하고,
2. SBP operator의 대각 질량 행렬과 경계 보간 특성을 이용해 각 요소에서 기울기를 인접 요소 경계값만으로 국소적으로 계산하는 방법을 제시하며,
3. 이를 통해 해와 기울기 모두에 대해 최적의 $p+1$ 차 수렴률(특히 기울기는 superconvergence)을 보이는 새로운 이산화 기법을 구현하였습니다.
이러한 접근법은 구현의 단순화뿐만 아니라, 고차 정확도를 요구하는 다양한 과학·공학 문제에 효과적으로 적용될 수 있습니다.