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

젖은 셀과 마른 셀을 한 시스템에 담고도 단조 수렴이 보장되는 이유.

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

논문 정보#

저자: 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이 자주 죽는가#

문제 형태는 단순하다.

\mathbf{V}(\eta) + T\eta = \mathbf{b}

여기서 $\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 분해).

a_i(\eta) = p_i(\eta) - q_i(\eta), \quad p_i, q_i \ge 0,\ \text{비감소}

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

V_i(\eta) = V_{1,i}(\eta) - V_{2,i}(\eta), \quad V_{j,i}(\eta) = \int_{-\infty}^\eta p_i\ \text{or}\ q_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. 안과 밖, 두 번의 선형화#

이제 다음 시스템을 푼다.

V_1(\eta) - V_2(\eta) + T\eta = \mathbf{b}

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

V_1(\eta^n) + (T - Q^{n-1})\eta^n = \mathbf{b} + V_2^{n-1} - Q^{n-1}\eta^{n-1}

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

(T + P^{n,m-1} - Q^{n-1})\eta^{n,m} = P^{n,m-1}\eta^{n,m-1} - V_1^{n,m-1} + \mathbf{d}^{n-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의 "초기치가 정답에 가까울 때만"이라는 약점이 사라진다. 아래 시뮬레이션에서 직접 조작해보자.

T = 0.50b = 0.90

step 1/4 · outer n=0 · eta=-0.4000 · |r|=1.10e+0

Pink = outer iterate (rises from below); green = inner iterate (falls from above). They squeeze toward the true root of V(eta)+T*eta=b with V(eta)=max(0,min(1,eta)).

T를 0.05까지 내리면 $T+P-Q$ 가 거의 0이 되어 정답이 갑자기 멀리 잡히고, 반복은 그래도 5번 안에 끝난다. b를 1.5 이상으로 올리면 정답이 $\eta > 1$ 구간(즉 $V$ 가 평탄해진 곳)에 떨어지면서 — Newton이라면 무한 진동할 자리에서 — 안쪽이 한 번에 멈춘다.

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

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

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

위 알고리듬을 1D 우물 단면에 그대로 돌려 매일 진화시키는 인터랙티브 데모다. 펌프 강도 슬라이더와 "+1 day" 버튼을 조작해본다.

pump flux phi = 0.0008 (m/s per cell)

day = 0 · pressurized cells = 58 · phreatic = 0 · dry = 2 · total water V = 3998.9 m

Cyan = pressurized (confined) region; teal = phreatic; black = dry. Yellow line is piezometric head eta(x). Pump near the center steadily drains the aquifer; the wet/dry boundary advances inward as days pass.

펌프를 0.0025 이상으로 올리고 10일 정도 돌리면 가운데에 검은 마른 구역이 나타나 점점 확장된다. 풀린 시스템은 매일 다른 활성 셀 수(active control volumes)를 가지지만 — Newton 한 번이라면 격자가 바뀌었다고 야코비안을 새로 잡아야 할 텐데 — 중첩 알고리듬은 같은 코드로 계속 풀린다.

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

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

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

재현 가능성 점수#

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

다음에 읽을 것#

Brugnano & Casulli (2008/2009) — piecewise linear systems의 Newton 유한 횟수 종료 정리. 이 논문의 직접 선조다.
Casulli (2009) — 같은 저자의 high-resolution wetting/drying 알고리듬.
Casulli & Zanolli (2010) — Richards 방정식 mixed form의 nested Newton 적용 (a(η)가 매끄러운 경우).