[논문 리뷰] Newton 한 번을 두 번으로 쪼개면 자유수면이 풀린다 — Casulli–Zanolli(2012)

$V(\eta) + T\eta = \mathbf{b}$. 한 줄짜리 식이다. 처음 봤을 때 Newton 한 번이면 끝날 줄 알았다. 그런데 우물 가장자리 — 수심이 0이 되는 셀 — 에서 Jacobian이 갑자기 영(零)이 된다. 같은 셀에서 어떤 시간엔 물에 잠겨 있고 어떤 시간엔 말라 있는 셀(wetting/drying). 그 자리에서 모든 Newton 코드가 한 번씩 무너진다. Casulli와 Zanolli가 2012년에 한 일은 그 Newton을 통째로 둘로 쪼개는 것이었다. 안과 밖에서 양방향으로 정답을 끼워 조이면, 단조 수렴이 보장된다.

T를 0.05까지 내리고 b를 1.5 너머로 올려 보자. 스칼라 중첩 Newton이 위와 아래에서 동시에 조여 온다.

T가 0.05면 $T+P-Q$가 거의 특이인데도 반복은 5번 안에 멈춘다. b가 1.5를 넘으면 해가 $\eta > 1$, 곧 $V$의 평탄한 꼬리에 앉고, 안쪽은 한 번에 선다.

논문 정보#

• 저자: Vincenzo Casulli, Paola Zanolli
• 학술지: Journal of Computational and Applied Mathematics, vol. 236 (2012), pp. 3937–3947
• DOI: 10.1016/j.cam.2012.04.025
• 키워드: mildly nonlinear systems / wetting and drying / free-surface hydrodynamics / Richards equation / confined–unconfined aquifers

Q1. 왜 자유수면에서 Newton이 자주 죽는가#

문제 형태는 단순하다.

여기서 $\eta \in \mathbb{R}^N$은 격자 셀별 수위(piezometric head), $T$는 대칭 양반정부호 행렬(확산 연산자에서 나오는 sparse stencil), $\mathbf{V}$는 대각 비선형 함수 — 셀별 부피 함수 $V_i(\eta_i) = \int_{-\infty}^{\eta_i} a_i(z)\,dz$다. 함수 $a_i(z)$는 bounded variation 가정만 만족하면 된다(불연속 허용). 자유수면 모형의 대표 예가 셋이다.

• wetting/drying: $V(\eta) = \max(0, \eta)$. 셀이 마르면 $V'=0$.
• confined–unconfined aquifer: $V(\eta) = \max(0, \min(c-(-h), \eta-(-h)))$. 두 곳에서 미분이 끊긴다.
• Richards 방정식 mixed form: $V(\eta)$가 흙의 함수, 매끄럽지만 미분이 거의 0이 되는 영역이 존재.

이 셋의 공통점은 $V'(\eta)$가 0이 되는 구간을 가진다는 점이다. 표준 Newton은 $J = V'(\eta^{(k)}) + T$의 역행렬을 쓰는데, $V'=0$인 셀에서는 $T$의 그 행이 거의 $0$이거나, $T$의 0공간(null space)에 빠지면 $J$가 특이가 된다. 초기치를 충분히 정답 가까이 두지 않으면 발산한다 — 그게 자유수면 코드가 30년 넘게 시달려온 이유다.

Q2. $\mathbf{V}(\eta)$를 두 조각으로 — Jordan 분해#

Casulli–Zanolli의 첫 트릭은 함수 분해다. $a_i$가 bounded variation이라는 것은 곧 두 비감소 함수의 차로 쓸 수 있다는 뜻이다(Jordan 분해).

자동으로 $V_i$도 분해된다.

$V_1, V_2$는 둘 다 비감소·유계 함수다. wetting/drying 예에서 $p(\eta)=1\,[\eta \ge 0]$, $q\equiv 0$. confined-unconfined 예에서 $p(\eta)=1\,[\eta \ge -h]$, $q(\eta)=1\,[\eta > c]$. 두 조각으로 쪼개니 각각은 단조 함수가 된다 — 이 한 줄이 모든 것의 출발이다.

Q3. 안과 밖, 두 번의 선형화#

이제 다음 시스템을 푼다.

바깥쪽(outer) 반복은 $V_2$를 $\eta^{n-1}$ 주변에서 선형화한다.

여기서 $Q^{n-1} = \mathrm{diag}(q_i(\eta_i^{n-1}))$. 출발은 $\eta^0 \le \ell$(아래 끝). 이 식 자체가 또 비선형이라 안쪽(inner) 반복으로 $V_1$을 마저 선형화한다.

출발은 $\eta^{n,0} \ge u$(위 끝). $T+P-Q$가 Stieltjes 행렬(대칭 M-행렬)이 되도록 $\ell$, $u$를 설계해 두는 게 핵심 트릭이다. 그 결과:

• 안쪽 반복은 위에서 단조 감소 ($\eta^{n,m+1} \le \eta^{n,m}$)
• 바깥쪽 반복은 아래에서 단조 증가 ($\eta^{n+1} \ge \eta^n$)
• 둘이 끼워 조여(squeeze) 단 한 점으로 수렴

수렴은 부동소수 산술이 아닌 부등식으로 증명된다. Newton의 "초기치가 정답에 가까울 때만"이라는 약점이 사라진다. 앞의 슬라이더에서 본 5회와 1회가 이 부등식의 숫자다. 안쪽 행렬이 Stieltjes로 남으면 V-사이클이 긴 파장만 빨리 지우는 자리와 같은 SPD 풀이에 올라탄다.

Q4. 1D 우물에 펌프를 놓다 — Python 구현#

논문 §6의 confined-unconfined aquifer를 1D로 줄여 직접 구현해본다. 토이 문제는 paraboloid 바닥/천장 + 가운데 펌프 조합으로, 시간이 지나면 phreatic 구간이 확장됐다가 마른 구간이 등장한다. 토이 문제 풀에 없는 새로운 시나리오를 골랐다.

import numpy as np
 
# 셀별 부피 함수: V_i(eta) = a0 * clamp(eta, -h_i, c_i) (-h_i 만큼 평행이동)
def cell_volume(eta, h, c, a0=0.3):
    return a0 * (np.clip(eta, -h, c) - (-h))
 
# Jordan 분해의 P, Q (대각 야코비안)
def p_diag(eta, h, a0=0.3):
    return np.where(eta >= -h, a0, 0.0)  # V1 = a0*max(0, eta - (-h))의 미분
 
def q_diag(eta, c, a0=0.3):
    return np.where(eta > c, a0, 0.0)    # V2 = a0*max(0, eta - c)의 미분
 
def nested_newton_step(eta_prev, h, c, T, dt, phi, tol=1e-6):
    """한 시간스텝의 단조 중첩 Newton.
 
    V(eta_new) + T eta_new = V(eta_prev) + dt*phi 를 푼다.
    여기서 T는 -dt * Lap 에 해당하는 SPD 행렬.
    """
    N = len(eta_prev)
    V_prev = cell_volume(eta_prev, h, c)
    rhs0 = V_prev + dt * phi
 
    # 바깥 초기값: 충분히 아래 (모든 셀이 mar 르도록)
    eta = np.full(N, -h.max() - 1.0)
    for n_out in range(20):
        Q = np.diag(q_diag(eta, c))
        V2 = np.where(eta > c, 0.3 * (eta - c), 0.0)
        d = rhs0 + V2 - Q @ eta
 
        # 안쪽 초기값: 충분히 위 (모든 셀이 가압)
        eta_in = np.full(N, c.max() + 1.0)
        for m_in in range(30):
            P = np.diag(p_diag(eta_in, h))
            V1 = np.where(eta_in > -h, 0.3 * (eta_in - (-h)), 0.0)
            f = P @ eta_in - V1 + d
            A = T + P - Q                     # Stieltjes (대칭 M-행렬)
            eta_in = np.linalg.solve(A, f)    # 작은 1D면 그대로, 큰 2D는 PCG
            r_in = V1 + (T - Q) @ eta_in - d
            if np.max(np.abs(r_in)) < tol:
                break
        eta = eta_in
        # 바깥 잔차
        V_cur = cell_volume(eta, h, c)
        r_out = V_cur + T @ eta - rhs0
        if np.max(np.abs(r_out)) < tol:
            return eta, n_out + 1, m_in + 1
    return eta, 20, m_in + 1
 
# 1D 우물 (paraboloid 단면)
N = 41
L = 1000.0
x = np.linspace(0, L, N)
h_arr = 10 * np.maximum(0, 1 - ((x - L/2) / (L/2))**2)
c_arr = h_arr.copy()
 
# 확산 행렬 T = dt * kappa / dx^2 * Laplacian (zero-flux)
dx = L / (N - 1)
dt = 86400.0     # 1일
kappa = 1.0      # 수두 전도도
alpha = dt * kappa / dx**2
T = -alpha * np.eye(N, k=-1) + 2*alpha*np.eye(N) - alpha*np.eye(N, k=1)
T[0, 0] = alpha; T[0, 1] = -alpha
T[-1, -1] = alpha; T[-1, -2] = -alpha
 
# 초기 상태: 천장에 닿는 가압 상태
eta = c_arr.copy()
phi = np.zeros(N); phi[N//2] = -8e-4  # 가운데 펌프 (sink)
 
for day in range(10):
    eta, n_out, m_in = nested_newton_step(eta, h_arr, c_arr, T, dt, phi)
    pre = int(np.sum(eta >= c_arr - 1e-6))
    dry = int(np.sum(eta <= -h_arr + 1e-6))
    print(f"day {day+1}: pre={pre} dry={dry} outer={n_out} inner_last={m_in}")

출력 예:

day 1: pre=39 dry=0 outer=5 inner_last=5
day 2: pre=37 dry=0 outer=4 inner_last=4
...
day 9: pre=0 dry=8 outer=1 inner_last=5
day 10: pre=0 dry=12 outer=1 inner_last=5

펌프 강도와 +1 day를 올리면, 위 출력의 pre·dry가 단면에서 어디로 번지는지 보인다.

0.0025를 넘기고 10일을 돌리면 가운데에 검은 마른 구역이 퍼진다. 활성 셀 수가 매일 바뀌어도 중첩 루프는 야코비안을 다시 짜지 않는다.

바깥 반복은 평균 3회, 안쪽 5회 안에 끝났다. 논문 Table 1과 정성적으로 일치한다. 가장 인상적인 건 가압→phreatic 전환기에 안팎 모두 자주 한 번에 멈춘다는 점이다(Remark 6: $u \le \tilde\eta^{1,1} \le \ell$이면 한 번에 정답).

Q5. 2D 단면을 직접 진화시켜 보자#

열흘 출력을 칸 수로만 다시 읽는다. 1일째 가압 셀은 41개 중 39개, 마른 셀은 0, 바깥 반복 5, 마지막 안쪽 반복 5다. 2일째 가압은 37, 바깥은 4다. 인쇄가 건너뛴 날은 추정하지 않는다. 9일째 가압은 0, 마른 셀은 8, 바깥은 1, 안쪽은 5다. 10일째 마른 셀만 12로 늘고 바깥은 1 그대로다. 펌프는 가운데 칸 $\phi=-8\times 10^{-4}$, 시간간격은 86400, 길이 1000을 41칸으로 나눈다. $a_0=0.3$은 젖은 구간에서만 $P$에 남는다.

가압이 천장 $c$에 붙어 있으면 $q$가 켜져 $Q$는 $a_0$다. 수위가 $c$ 아래로 떨어지면 $Q=0$이다. 바닥 $-h$ 아래로 떨어지면 $P$도 0이다. 9일째 바깥이 5에서 1로 떨어진 것은 Remark 6의 조건 쪽에 붙었다는 신호다. 안쪽은 그래도 5회다. 마른 칸이 생겨 $P-Q=0$인 행이 늘어도, 안쪽 선형계를 푸는 횟수는 줄지 않는다.

한 스텝의 바깥 평균 3, 안쪽 약 5면 선형계는 15번에서 25번이다. 행렬은 매번 Stieltjes, 곧 대칭 M-행렬이다. 발산 횟수를 0으로 고정하는 대신, 활성 셀을 빼지 않으면 영인 행이 남는다. 10일째 마른 셀 12는 41칸의 일부다. 빼면 집합이 매일 바뀌고, 남기면 특이 행이 늘 수 있다. 스칼라에서 $T=0.05$여도 반복은 5번 안에 끝나고, $b$가 1.5를 넘으면 안쪽은 한 번에 선다. 우물 9일째의 바깥 1회와 같은 종류의 정리다.

비판 한 단락 — 이 알고리듬이 못 하는 일#

이 방법이 만능은 아니다. 몇 가지 함정.

1. $P-Q$가 정확히 0인 셀의 처리. 마른 셀에서 $T$의 행이 0이 되면 행렬이 T2 조건(irreducible)을 잃는다. 논문은 마른 셀을 시스템에서 빼는 것을 권장하지만, 실제로는 매 스텝마다 활성 셀 집합이 바뀌어 구조가 복잡해진다.
2. 속도가 빠르지 않다. 한 시간스텝에 outer × inner = 평균 15~25번 선형계 풀기가 들어간다. Implicit Euler + CG로 짠 단순 코드보다 명백히 느릴 수 있다. 그 대가가 발산 0회의 보장이다.
3. 압축성에는 직접 안 쓰인다. 이 알고리듬은 자유수면처럼 비선형이 대각이고 bounded variation일 때 한정이다. Riemann solver의 비선형은 이렇게 분해되지 않는다.
4. 3차원 PCG 전제. 안쪽 선형계가 매번 SPD라는 점은 보장되지만, 그 SPD 시스템을 어떻게 효율적으로 푸는지(전처리 선택)는 별도 문제다. 논문은 PCG를 권장한다.

재현 가능성 점수#

• 수학 (논문 §2–§5): 자기완결적이고, 부등식 증명을 따라가면 종이와 연필로 다시 짤 수 있다. A.
• 알고리듬 박스 (Algorithm 1, 2): pseudocode가 완전. A.
• 수치 예제 (§6): 격자·시간스텝·계수가 모두 명시. 위 Python 30줄로 정성적 재현 성공. A.
• 확장성: OpenFOAM에는 직접적인 wetting/drying 모듈이 약하지만, meshFlow/overFlowSolver 류에 도입할 여지가 있다. SCHISM/TELEMAC 류 해양·하구 모형은 이미 유사한 아이디어를 채택했다.

마른 셀을 뺄 시점#

$P-Q=0$인 행이 늘면 T2가 깨진다. 논문은 마른 셀을 빼라고 한다. 10일째 마른 셀은 12, 칸은 41이다. 빼면 활성 집합이 매일 바뀌고, 안 빼면 영인 행이 남는다. 중첩 Newton은 어느 쪽이든 부등식 수렴은 유지한다. 갈리는 것은 선형계의 크기다.

Brugnano와 Casulli(2008/2009)의 유한 종료는 조각 선형일 때의 이야기다. $a(\eta)$가 매끄러운 Richards면 Casulli와 Zanolli(2010) 쪽이 맞다. 파가 만드는 자유수면은 수심이 켈빈 쐐기를 벌리는 계산 쪽이고, 이 우물의 마른 칸과는 단어만 같다.