불연속 갤러킨법

편집 역사 토론
수치해석 전산유체역학 마지막 수정: 2026-08-03 04:51:44

1. 개요[편집]

유한요소법의 고차 다항식과 유한체적법의 리만 솔버를 한 그릇에 넣고 비볐더니, 질량행렬이 블록 대각으로 떨어졌다.

불연속 갤러킨법(Discontinuous Galerkin method, DG)은 각 요소 내부에서 해를 다항식으로 전개하되 요소 경계에서 연속성을 전혀 강제하지 않고, 인접 요소와의 결합을 오직 경계면의 수치 플럭스로만 처리하는 이산화 기법이다. 갤러킨 방법의 약형식 골격은 그대로 쓰지만, 시험함수 공간을 전역 연속 공간이 아닌 요소별 독립 다항식 공간 Vh={v:v∣K∈Pp(K)}V_h = \{v : v|_K \in \mathbb{P}^p(K)\} 으로 잡는다는 점이 다르다. 요소 하나가 곧 독립된 작은 유한요소법 문제이고, 그 사이를 유한체적법의 방식으로 잇는 셈이다.

1973년 리드(Reed)와 힐(Hill)이 중성자 수송 방정식용으로 처음 제안했고, 1989~1998년 콕번(Cockburn)과 슈(Shu)의 RKDG 연속 논문으로 비선형 보존 법칙에 대한 현대적 형태가 확립됐다.1 오늘날 고차 압축성 유동, 전산음향학, 대와류 모사, 전자기 시간영역 해석에서 주력으로 쓰인다.

2. 약형식과 수치 플럭스[편집]

보존형 방정식 ∂tu+∇⋅f(u)=0\partial_t u + \nabla \cdot \mathbf{f}(u) = 0 에 요소 KK 국소로 시험함수 ϕ\phi 를 곱하고 부분적분하면

∫K∂uh∂t ϕ dV  −  ∫Kf(uh)⋅∇ϕ dV  +  ∮∂Kf^(uh−,uh+)⋅n ϕ dS=0\int_K \frac{\partial u_h}{\partial t}\,\phi \,dV \;-\; \int_K \mathbf{f}(u_h)\cdot\nabla\phi \,dV \;+\; \oint_{\partial K} \hat{\mathbf{f}}(u_h^-, u_h^+)\cdot \mathbf{n}\,\phi \,dS = 0

가 된다. 여기서 uh−u_h^-, uh+u_h^+ 는 경계면 양쪽에서 본 값이고 일반적으로 서로 다르다. 그 불일치를 하나의 값으로 정리하는 함수 f^\hat{\mathbf{f}} 가 수치 플럭스이며, 이 자리에 유한체적법에서 쓰던 리만 솔버를 그대로 꽂아 넣으면 된다. 업윈드, 로컬 락스-프리드리히스(Rusanov), HLL/HLLC, 로(Roe) 무엇이든 좋다.

f^LLF=12(f(u−)+f(u+))⋅n−12α (u+−u−),α=max⁡∣λ∣\hat{\mathbf{f}}^{\text{LLF}} = \tfrac{1}{2}\big(\mathbf{f}(u^-) + \mathbf{f}(u^+)\big)\cdot\mathbf{n} - \tfrac{1}{2}\alpha\,(u^+ - u^-), \qquad \alpha = \max|\lambda|

이 재활용성이 DG의 최대 정치적 자산이다. FVM 커뮤니티가 30년간 쌓은 근사 리만 솔버 자산을 하나도 안 버리고 차수만 올릴 수 있다. 플럭스가 일치성(consistency)과 단조성/E-flux 조건을 만족하면 국소 보존은 요소 단위로 정확히 유지된다.

3. 위치 — FVM과 FEM 사이[편집]

항목FVMDG연속 FEM
요소 간 연속성없음(셀 평균)없음C0C^0 강제
고차화 방법이웃 셀 재구성(스텐실 확장)요소 내부 자유도 추가차수 상승
병렬 통신넓은 헤일로면 이웃 1겹절점 공유
질량행렬대각(자명)블록 대각전역 희소
자유도 수셀당 1셀당 (p+1)d(p{+}1)^d절점 공유로 절약

FVM에서 3차 이상으로 가려면 MUSCL/WENO 재구성으로 스텐실을 이웃의 이웃까지 넓혀야 하고, 비정렬 격자에서는 이게 지옥이다. DG는 차수를 올려도 스텐실이 면 이웃 1겹으로 고정된다. 통신량이 차수와 무관하니 대규모 병렬에서 스케일링이 좋고, GPU와도 궁합이 좋다(GPU 컴퓨팅).

반대로 연속 FEM과 비교하면 절점을 공유하지 않아 자유도가 늘어난다. 3차원 p=3p=3 육면체 기준 DG는 요소당 64개 자유도를 홀로 들고 있고, 연속 FEM은 이웃과 나눠 쓴다. 저차에서는 DG가 명백히 비싸고, 고차로 갈수록 그 비율이 좁혀지며 정확도로 역전한다. 그래서 DG는 p≥2p \ge 2 부터가 사는 영역이다.

4. 블록 대각 질량행렬 — 명시적 적분의 공짜 점심[편집]

DG의 자유도는 요소 밖으로 나가지 않으므로 질량행렬

MijK=∫Kϕiϕj dVM^K_{ij} = \int_K \phi_i \phi_j \,dV

이 요소 단위 블록 대각이다. 전역 연립방정식을 푸는 게 아니라 (p+1)d(p{+}1)^d 크기의 작은 블록을 각자 역행렬 하면 끝이고, 그마저도 직교다항식(르장드르 계열) 모달 기저함수를 쓰고 아핀 사상 요소면 대각으로 떨어진다. 즉 명시적 시간적분에서 질량행렬 역이 사실상 공짜다.

이 성질 때문에 DG는 SSP(strong stability preserving) 룽게-쿠타법과 결합해 완전 명시적 스킴으로 굴리는 것이 표준이 됐다. 연속 FEM은 같은 자리에서 전역 질량행렬을 만나 mass lumping이라는 정확도 손실을 감수해야 한다.

또 하나의 자연스러운 이점이 p-적응이다. 요소마다 다항식 차수를 독립적으로 바꿔도 인접 요소와의 연속성을 맞출 필요가 없다 — 어차피 안 맞추니까. hanging node도 마찬가지로 자유로워서 적응 격자 세분화와 결합한 hp-적응이 구현상 거의 공짜다. 매끄러운 해에서 hp-적응은 지수 수렴을, 특이점 근처에서는 h-세분으로 대수 수렴을 준다.

5. 대가 — CFL과 진동[편집]

공짜는 없다. 명시적 DG의 시간 스텝 제한은

Δt  ≤  C h(2p+1) ∣λ∣max⁡\Delta t \;\le\; C\,\frac{h}{(2p+1)\,|\lambda|_{\max}}

로, 차수 pp 에 대해 1/(2p+1)1/(2p+1) 로 조여진다.2 p=4p=4 면 같은 격자에서 CFL 조건이 1차 스킴의 1/9이다. 게다가 pp 를 올리면 자유도도 늘어나므로 벽시계 시간은 이중으로 나빠진다. DG가 “정확도당 비용”에서 이기는 것이지 “격자당 비용”에서 이기는 게 아니라는 점은 늘 헷갈리는 대목이다. 이 제약을 피하려고 요소별로 다른 시간 스텝을 쓰는 국소 시간적분(LTS)이나 암시적/IMEX 적분을 붙이는데, 암시적으로 가는 순간 블록 대각의 축복은 사라진다.

두 번째 대가는 불연속 근처의 진동이다. 고차 다항식으로 충격파를 표현하면 깁스 현상이 그대로 나와 밀도가 음수가 되고 계산이 죽는다. 대응책은 하나같이 불완전하다.3

  • 제한자(slope/moment limiter): minmod 기반 TVB 제한자(차분 도식의 TVD 논리와 같은 뿌리), 모멘트 제한자, WENO 도식 기반 제한자. 문제는 매끄러운 극값에서도 발동해 차수를 1차로 떨어뜨리는 것(clipping). TVB 상수 MM 을 사용자가 튜닝해야 하는데, 이 값이 문제 의존적이라 국룰이 없다.
  • 인공점성: 페르손-페라이레(2006)의 매끄러움 지시자로 문제 요소를 골라 요소별 점성 ε∼h/p\varepsilon \sim h/p 를 넣는다. 충격을 요소 폭 안에 뭉개는 방식이라 매끄럽지만, 계수 튜닝은 여전히 남는다.
  • 하이브리드: 문제 요소만 1차 FV 서브셀로 갈아 끼우는 서브셀 유한체적 제한(subcell FV limiting). 요즘 압축성 코드에서 선호도가 높다.
  • 양수성 보존 제한자: 장-슈(Zhang–Shu, 2010)의 스케일링 제한자는 셀 평균이 유효하면 밀도·압력의 양수성을 증명 가능하게 보존한다. 진동을 없애진 않지만 코드가 죽지는 않게 해 준다.

여기에 비선형 항의 앨리어싱도 있다. 부족 적분(under-integration)으로 생기는 앨리어싱 오차는 고차에서 발산을 유발하며, 과적분(over-integration)이나 분할형(split-form)·엔트로피 안정 정식화로 대응한다.

6. 타원 문제와 IPDG[편집]

이류 문제는 플럭스로 잘 붙지만, 라플라스 방정식 같은 2계 타원 문제는 그냥 두면 요소 간 점프를 아무도 벌하지 않아 문제가 잘 정의되지 않는다. 그래서 내부 벌점법(Interior Penalty DG, IPDG)이 쓰인다. 대칭 IPDG(SIPG)의 쌍선형 형식은 요소 적분에 더해 각 면에서

−∫e{∇uh}⋅n [vh]  −  ∫e{∇vh}⋅n [uh]  +  σhe∫e[uh][vh]- \int_e \{\nabla u_h\}\cdot\mathbf{n}\,[v_h] \; - \; \int_e \{\nabla v_h\}\cdot\mathbf{n}\,[u_h] \; + \; \frac{\sigma}{h_e}\int_e [u_h][v_h]

세 항을 더한다. 앞 두 항이 일치성과 대칭성을 주고, 마지막 벌점항이 점프를 억제한다. 여기서 벌점 상수 σ\sigma 가 생명줄이다. 너무 작으면 강제성(coercivity)을 잃어 시스템이 부정부호가 되고, 너무 크면 조건수가 σp2/h\sigma p^2/h 로 치솟아 반복 solver가 죽는다. 이론적 하한이 σ≳Cp2\sigma \gtrsim C p^2 로 알려져 있고, 실무에서는 그 근처의 최소값을 쓰는 것이 정석이다.4 벌점항 없이도 되는 변형으로 LDG, BR2 등이 있으며, 압축성 점성항에서는 BR2 계열이 널리 쓰인다.

7. 관련 문서[편집]

8. Footnotes[편집]

  1. Reed, W. H. & Hill, T. R. (1973), LA-UR-73-479. 중성자 수송용으로 만든 방법이 20년 뒤 CFD의 주력이 된 셈인데, 원래 문제가 선형 이류라서 수치 플럭스가 그냥 업윈드였다는 게 포인트. 비선형 보존 법칙으로 넘어오면서 어려운 부분이 전부 플럭스와 제한자로 몰렸다. ↩

  2. Cockburn & Shu의 RKDG 해석에서 나온 결과로, 1/(2p+1)1/(2p+1) 은 1차원 선형 이류 + SSP-RK 기준의 값이다. 다차원·비정렬 격자에서는 요소 형상과 최소 내접구 반경이 개입해 상수가 더 나빠진다. 요컨대 코드에 박힌 CFL 수는 항상 경험적으로 다시 잡아야 한다. ↩

  3. 그래서 고차 코드 튜토리얼의 절반이 “왜 내 계산이 3000스텝에서 NaN이 나나요”에 대한 답이다. 답은 대체로 제한자를 안 켰거나, 켰는데 TVB 상수가 문제에 안 맞거나, 적분 차수가 모자라서 앨리어싱이 터진 것 셋 중 하나. ↩

  4. 벌점 상수를 “일단 1000쯤 박으면 안전하겠지”라고 하는 순간 반복 solver의 수렴이 죽는다. 강제성은 확보되는데 조건수가 같이 올라가기 때문. 안정성과 조건수를 동시에 만족시키는 최소 σ\sigma 를 찾는 것이 IPDG 튜닝의 전부라고 해도 과언이 아니다. ↩