MUSCL

1. 개요

MUSCL Monotone Upstream-centred Schemes for Conservation Laws
고안	Bram van Leer (1973~1979)
핵심	셀 내부를 조각선형으로 재구성 → 공간 2차
진동 억제	기울기 제한자(minmod·van Leer·superbee)로 TVD 유지
틀	재구성 → 발전 → 평균 (REA)
기반	고두노프 도식 + 리만 솔버

셀 안이 상수라고 우기던 고두노프에게 “기울기 하나만 얹자”고 제안한 게 이 바닥을 2차 시대로 밀어 넣었다.

MUSCL(Monotone Upstream-centred Schemes for Conservation Laws)¹²은 유한체적법에서 각 셀의 값을 상수가 아니라 조각선형(piecewise-linear)으로 재구성해 공간 정확도를 1차에서 2차로 올리되, 재구성 기울기에 제한자(limiter)를 물려 불연속 근처의 진동을 막는 고해상도 도식(高解像度, high-resolution scheme)의 원형이다. 네덜란드 출신 수치해석가 브람 판 레이어(Bram van Leer)가 1970년대에 발표한 다섯 편짜리 연작 Towards the Ultimate Conservative Difference Scheme I–V(1973–1979)에서 골격이 잡혔고, 특히 마지막 5편(1979)이 오늘날 MUSCL이라 불리는 그 도식이다.³

핵심 발상은 한 문장이다. 고두노프 도식이 셀을 계단으로 보는 대신, 셀마다 기울어진 직선으로 보라. 계단은 1차, 직선은 2차다. 그런데 직선을 아무렇게나 그으면 급구배 근처에서 오버슈트가 터지므로, 기울기를 해에 따라 눌러 주는 비선형 장치가 반드시 따라붙는다. 그 장치가 제한자이고, 왜 하필 비선형이어야 하는가(고두노프 장벽)와 제한자 설계 규격(스위비 다이어그램)은 전변분 감소 문서가 소스 오브 트루스다 — 여기서는 반복하지 않고 재구성이라는 MUSCL 고유의 각도만 판다.

2. 왜 재구성인가 — 고두노프의 한계

리만 솔버 문서에서 보듯, 고두노프 법은 셀 값을 상수로 두어 모든 셀 경계를 리만 문제로 만든다. 충격파를 물리적으로 옳게 잡는다는 큰 장점이 있지만, 셀 안을 상수로 뭉갠 대가로 공간 1차에 머문다. 1차 상류 도식의 병폐는 수치 소산과 분산에서 다루는 그대로다 — 접촉면과 매끈한 경사가 몇 셀에 걸쳐 번져 버린다.

정확도를 올리는 정직한 길은 경계면에 넘겨줄 좌·우 상태값을 셀 평균 하나가 아니라 이웃 정보로 보간해 더 잘 추정하는 것이다. 셀 $i$ 안의 해를

$$u_i(x) = \bar{u}_i + \sigma_i\,(x - x_i), \qquad x \in [x_{i-1/2},\, x_{i+1/2}]$$

인 직선으로 재구성하면, 셀 평균 $\bar u_i$ 는 그대로 보존($\int$ 이 셀 중심에서 상쇄)하면서 경계에서의 값이

$$u_{i}^{R} = \bar u_i + \tfrac{1}{2}\sigma_i \Delta x, \qquad u_{i}^{L} = \bar u_i - \tfrac{1}{2}\sigma_i \Delta x$$

로 갱신된다. 기울기 $\sigma_i$ 를 중심차분 $(\bar u_{i+1}-\bar u_{i-1})/2\Delta x$ 로 잡으면 매끈한 곳에서 2차다. 문제는 이 중심 기울기가 불연속 근처에서 언더슈트·오버슈트를 낳는다는 것 — 그래서 기울기를 그대로 쓰지 않고 제한한다.

3. 기울기 제한 — 플럭스 제한자와의 쌍대

MUSCL은 제한자를 기울기 자체에 건다(slope limiter). 상류·하류 차분

$$a = \frac{\bar u_i - \bar u_{i-1}}{\Delta x}, \qquad b = \frac{\bar u_{i+1} - \bar u_i}{\Delta x}$$

두 개를 놓고 제한된 기울기를 예컨대 minmod로

$$\sigma_i = \operatorname{minmod}(a,\,b) = \begin{cases} a & |a|<|b|,\ ab>0 \\ b & |b|\le|a|,\ ab>0 \\ 0 & ab\le 0 \end{cases}$$

로 정한다. $ab\le 0$ 은 셀 $i$ 가 극값이라는 뜻이고, 그때 $\sigma_i=0$ 이면 그 셀은 1차 상류로 내려앉는다 — 극값에서 1차로 후퇴하는 이 동작이 진동을 막는 대가이자, 전변분 감소에서 말한 클리핑(극값 뭉개기)의 근원이다. minmod 대신 van Leer·van Albada·MC·superbee를 쓰면 같은 틀에서 소산의 강도만 달라진다.

여기서 짚어 둘 것: 선형 이류에 한해 기울기 제한자(MUSCL)와 플럭스 제한자(TVD 플럭스형)는 서로 변환되는 같은 것이다. 스위비 변수 $r$ 로 쓴 플럭스 제한자 $\psi(r)$ 과 MUSCL 기울기 비 $\sigma$ 사이에 일대일 대응이 있어, 같은 제한자가 두 언어로 나타날 뿐이다. 비선형 계나 다차원으로 가면 두 형식이 갈라지지만, 입문 단계에서 “MUSCL이냐 TVD냐”를 대립시키는 것은 대개 헛다리다.

4. 재구성-발전-평균 (REA)

MUSCL을 한 스텝 돌리는 절차를 르베크(LeVeque)는 재구성–발전–평균(Reconstruct–Evolve–Average, REA)이라는 세 박자로 정리했다.

1. 재구성(Reconstruct). 셀 평균 $\{\bar u_i\}$ 에서 제한된 조각선형 $u_i(x)$ 를 만든다. 여기서 경계 양쪽의 외삽값 $u_i^{R}$, $u_{i+1}^{L}$ 이 나온다.
2. 발전(Evolve). 각 경계에서 좌·우 외삽값을 초기조건으로 하는 리만 문제를 풀어 시간 전진시킨다. 실무에서는 정확 발전 대신 리만 솔버로 경계 플럭스 $F_{i+1/2}$ 만 뽑는다.
3. 평균(Average). 발전된 해를 다시 셀 평균으로 되돌려 $\bar u_i^{n+1}$ 을 얻는다. 이 단계가 보존형 갱신

$$\bar u_i^{n+1} = \bar u_i^{n} - \frac{\Delta t}{\Delta x}\bigl(F_{i+1/2} - F_{i-1/2}\bigr)$$

그 자체다.

이 틀의 값어치는 “고차 = 재구성을 더 정교하게”라는 한 축과 “물리 = 리만 솔버가 담당”이라는 다른 축을 깔끔히 분리한다는 데 있다. 재구성을 조각선형에서 조각포물선으로 올리면 조각포물선법(PPM, Colella–Woodward 1984)이 되고, 매끄러움 가중으로 갈아타면 WENO 도식과 ENO 도식이 된다. MUSCL은 이 사다리의 첫 칸이다.

5. MUSCL-Hancock — 반스텝 예측자

경계 상태를 시간으로도 2차로 만들려면 재구성한 값을 반 스텝 앞으로 밀어 놓고 리만 문제를 풀어야 한다. 판 레이어와 행콕(Hancock)의 방식이 표준이다.

$$\bar u_i^{L/R,\,*} = u_i^{L/R} - \frac{\Delta t}{2\,\Delta x}\Bigl(F(u_i^{R}) - F(u_i^{L})\Bigr)$$

즉 셀 내부 플럭스 차이로 좌·우 외삽값을 반 스텝만큼 자체 발전시킨 뒤(예측자, predictor), 이 예측값 $u_i^{R,*}$ 과 $u_{i+1}^{L,*}$ 을 경계 리만 솔버에 넣어 최종 플럭스를 얻는다(수정자, corrector). 리만 솔버를 한 스텝에 한 번만 호출하고도 시공간 2차를 얻는다는 점에서 계산 효율이 좋아, 압축성 코드에서 오래 사랑받았다.

시간 정확도를 얻는 다른 길은 공간만 MUSCL로 재구성해 반이산(semi-discrete) 형태로 두고, 시간을 SSP 룽게-쿠타(룽게-쿠타법의 강안정성 보존 변형)로 올리는 것이다. 2차·3차 SSP-RK가 CFL 계수 1로 나와 있어 사실상 공짜이고, 불연속 갤러킨법과 코드를 공유하기 좋아 요즘은 이쪽이 더 흔하다.

버거스 방정식 $u_t + (u^2/2)_x = 0$ 같은 스칼라 문제에서 초기 매끈한 프로파일이 자발적으로 충격을 만들 때, MUSCL은 충격을 CFL 조건 한도 안에서 두어 셀 폭으로 붙잡으면서 그 앞뒤의 매끈한 영역은 2차로 유지한다. 다만 충격 자체의 폭에서는 차이가 생각보다 작다 — 버거스의 볼록 플럭스가 특성선을 스스로 모으기 때문에 1차 고두노프도 충격을 3셀 안에 잡고 MUSCL은 2셀로 한 칸 줄일 뿐이다(CFL 0.20.8에서 동일). 재구성 차수가 진짜로 값을 하는 곳은 희박팬과 매끈한 극값이며, 위 시뮬의 희박 프리셋에서 $L^1$ 오차가 $1.13\times10^{-2}$에서 $1.4\text{–}2.7\times10^{-3}$로 68배 줄어드는 것이 그 증거다.

6. 검산 — 차수와 TVD를 직접 재 본다

MUSCL을 두고 “2차 정확도이면서 진동이 없다”고 뭉뚱그리는 설명이 많지만, 두 주장은 같은 제한자에서 동시에 성립하지 않는다. 선형 이류 $u_t + u_x = 0$을 주기 영역에서 MUSCL-행콕 한 스텝 형태로 풀어 직접 재면 바로 드러난다.

6.1. 매끈한 해에서의 관측 차수

초기조건 $u_0 = \sin 2\pi x$, CFL 0.4, $t=1$(한 주기)에서 $L^1$ 오차와 격자를 두 배씩 줄여 얻은 관측 차수.

도식	N=40	N=80	N=160	N=320	N=640	관측 차수
1차 풍상	1.63e−1	8.76e−2	4.54e−2	2.31e−2	1.17e−2	0.90 → 0.99
minmod	2.93e−2	8.32e−3	2.26e−3	6.19e−4	1.64e−4	1.81 → 1.91
van Leer	1.15e−2	2.78e−3	6.46e−4	1.47e−4	3.32e−5	2.05 → 2.15
MC	6.48e−3	1.51e−3	3.25e−4	6.81e−5	1.47e−5	2.11 → 2.21
superbee	1.78e−2	6.01e−3	1.69e−3	4.44e−4	1.13e−4	1.56 → 1.97
제한자 없음	1.38e−2	3.45e−3	8.64e−4	2.16e−4	5.40e−5	2.00 정확히

읽을 것이 둘이다. 첫째, 제한자를 끈 쪽만 차수가 깔끔하게 2.00이다. 제한자를 물린 순간 차수는 2에서 흐트러진다 — minmod는 1.9 근처에 머물고(사인파의 극값 두 곳에서 $ab \le 0$이 되어 $\sigma=0$으로 주저앉는, 바로 그 클리핑의 대가다), superbee는 거친 격자에서 1.56까지 내려간다. 즉 “TVD 제한자를 달면 2차”가 아니라 **“매끈한 단조 구간에서만 2차”**가 정확한 진술이다.⁴ 둘째, 그래도 제한자를 단 쪽이 1차 풍상보다 두 자릿수 가까이 정확하다 — $N=640$에서 minmod의 오차는 1차의 1/71이다.

6.2. 기울기 제한자 = 플럭스 제한자 (수치로 확인)

위에서 “선형 이류에 한해 두 형식은 같은 것”이라고 했는데, 이것은 비유가 아니라 대수적 동일성이다. 기울기 $\sigma_i = \phi(r_i)(u_{i+1}-u_i)$로 재구성한 뒤 반 스텝 발전시킨 좌측값은

$$u_i^{R,*} = u_i + \tfrac{1}{2}(1-\nu)\,\phi(r_i)\,(u_{i+1}-u_i), \qquad \nu = \frac{a\Delta t}{\Delta x}$$

이고, 이것이 곧 스위비의 플럭스 제한자 형태 $F_{i+1/2} = a u_i + \tfrac12 a(1-\nu)\phi(r_i)(u_{i+1}-u_i)$다. 두 코드를 따로 구현해 200셀을 한 주기 돌리면 차이가 배정도 반올림 수준에 머문다.

제한자	두 형식의 max|차이| (t=1, 200셀)
minmod	7.8×10⁻¹⁶
van Leer	6.7×10⁻¹⁶
superbee	8.9×10⁻¹⁶
MC	8.9×10⁻¹⁶

6.3. 단조성은 지켜지는가

사각파 $u_0 = 1_{[0.3,0.6)}$를 200셀에서 한 주기 이류시킨 뒤의 극값과 전변분. 정확해의 전변분은 2다.

도식	min u	max u	TV	L¹ 오차
minmod	0.0000	1.0000	2.0000	3.35e−2
van Leer	0.0000	1.0000	2.0000	2.14e−2
MC	0.0000	1.0000	2.0000	1.78e−2
superbee	0.0000	1.0000	2.0000	8.97e−3
제한자 없음	−0.2388	1.2386	4.2367	5.72e−2
1차 풍상	0.0000	0.9939	1.9877	8.74e−2

네 제한자 모두 전변분을 2.0000으로 정확히 유지하고 오버슈트가 0이다. 제한자를 끄면 전변분이 2.1배로 불어나고 24%짜리 언더슈트가 생긴다 — 위 차수 표에서 2.00을 기록한 그 도식이다. 차수와 단조성의 교환이 이 두 표 사이에 그대로 박혀 있는 셈이다. 사각파에서는 압축적인 superbee가 L¹ 오차를 minmod의 1/3.7로 줄이는데, 같은 superbee가 사인파에서는 가장 나쁜 차수를 보였다는 것도 함께 읽어야 한다. 제한자에 공짜 우승자는 없다.

참고로 2차 정확도의 필요조건은 $\phi(1) = 1$이고(매끈한 구간 $r \approx 1$에서 중심 기울기를 복원해야 한다), 네 제한자 모두 이를 정확히 만족한다. 스위비 다이어그램의 2차 TVD 영역이 minmod와 superbee 사이의 띠로 그려지는 것도 이 조건 때문이다.

r	minmod	van Leer	MC	superbee
0.25	0.250	0.400	0.500	0.500
0.5	0.500	0.667	0.750	1.000
1	1.000	1.000	1.000	1.000
2	1.000	1.333	1.500	2.000
10	1.000	1.818	2.000	2.000

7. 계와 다차원으로

스칼라 이류를 넘어 압축성 유동의 오일러 방정식 같은 계(system)로 가면 손봐야 할 곳이 생긴다. 보존 변수(밀도·운동량·에너지)에 성분별로 제한자를 그냥 먹이면 접촉면 근처에서 압력 진동이 새기 쉽다. 정석은 리만 솔버의 파동 분해에서 나온 특성 변수로 갈아탄 뒤 성분별로 재구성하는 것이다. 국소적으로 스칼라 문제 여러 개를 푸는 셈이라 진동이 잘 잡히지만, 매 셀 경계에서 좌·우 고유벡터를 곱하는 비용이 붙는다. 싸게 가려고 원시 변수(밀도·속도·압력)에 제한자를 먹이는 절충이 실무에서 흔하다.

다차원에서 진짜 2차 TVD는 사실상 불가능하다는 것도 전변분 감소에 적힌 그대로여서, MUSCL도 대개 방향 분할이나 차원별 재구성으로 우회한다. 강한 충격이 격자에 비스듬히 지날 때 방향 분할의 오차가 도드라지면, 진짜 다차원 재구성이나 적응 격자 세분화로 넘어간다.

8. 실무에서 남는 것

• 제한자 선택이 곧 소산-선명도 다이얼이다. 충격이 주인공이면 압축적인 제한자(MC·superbee 쪽), 난류·음향처럼 매끈한 구조가 중요하면 소산이 큰 minmod·van Leer 쪽. superbee는 매끈한 프로파일을 사각파로 만들어 놓으니 계면 추적이 아니면 피한다.
• 정상상태 수렴을 방해하는 제한자 채터링은 여기서도 나온다. 미분가능한 제한자나 벤카타크리슈난식 완화가 그 처방.
• 차수 검증은 매끈한 해로, 충격 포착은 불연속으로 따로 본다. 하나의 문제로 둘 다 재려 하면 클리핑에 차수가 잡아먹혀 진단이 꼬인다. 검증 및 확인의 기본기.

9. 관련 문서

• 유한체적법 · 리만 솔버 · 고두노프 도식
• 전변분 감소 · WENO 도식 · ENO 도식 · 조각포물선법
• 수치 소산과 분산 · 차분 도식 · CFL 조건
• 룽게-쿠타법 · 불연속 갤러킨법 · 적응 격자 세분화
• 압축성 유동 · 충격파 · 나비에-스토크스 방정식 · 전산유체역학

10. Footnotes

1. MUSCL의 M이 monotone인데 정작 만들어진 도식은 엄밀한 의미의 단조(monotone) 도식이 아니다 — 단조 도식은 1차가 한계이기 때문이다. 여기서 monotone은 “단조성을 보존하려 애쓰는” 정도의 느슨한 뜻이고, 정확히는 TVD다. 약칭은 대개 이렇게 슬쩍 과장된다. ↩

2. “upstream-centred”라는 말도 재밌다. 상류(upwind)로 정보를 가져오되 중심(centred) 재구성을 얹는다는, 두 진영을 다 먹겠다는 선언이다. 실제로 MUSCL의 성공은 상류의 안정성과 중심차분의 정확도를 한 도식에 절충해 넣은 데 있다. ↩

3. 연작의 제목 Towards the Ultimate Conservative Difference Scheme(“궁극의 보존형 차분 도식을 향하여”)은 수치해석 논문 제목 중 손꼽히는 허세이자, 실제로 그 목표를 거의 달성해 버린 드문 사례다. 판 레이어는 훗날 이 연작으로 이 분야의 고전 반열에 올랐고, MUSCL이라는 약칭도 5편에서 본인이 붙였다. ↩

4. 그래서 제한자를 단 도식의 차수 검증에는 함정이 있다. 극값이 없는 단조 프로파일(예: 한 방향으로만 증가하는 매끈한 층)로 재면 minmod도 2.0을 준다. 논문의 차수 표가 2.0으로 깔끔할 때 어떤 초기조건을 썼는지 확인해 볼 값어치가 여기 있다. ↩