[논문 리뷰] 같은 점의 색함수가 경로에 따라 0.544 달랐다 — 컬 구속조건과 엇갈린 격자

같은 점에서 색함수가 두 값으로 읽혔다#

계면 장 $\mathbf{b}$를 400스텝 굴린 뒤, 격자 한가운데 셀의 색함수를 복원해봤다. $\mathbf{b}$를 원점부터 적분하면 색함수가 나와야 한다. 정의가 $\mathbf{b} = \nabla c$이니까. $x$로 먼저 가고 $y$로 꺾은 값은 한 숫자였다. $y$로 먼저 가고 $x$로 꺾은 값은 다른 숫자였다. 차이가 0.544였다. 색함수는 0과 1 사이에 있는 양이다.

전단 세기 슬라이더와 Rusanov 소산 슬라이더를 같이 밀면, 동그란 계면이 한쪽 격자에서만 색을 띠기 시작한다.

왼쪽 패널이 엇갈린 배치, 오른쪽이 충돌 배치다. 왼쪽이 검게 머무는 동안 오른쪽이 언제 밝아지는지를 보면, 경로 차 0.544가 어느 쪽에 쌓이는지 보인다.

이 글은 Chiocchetti와 Dumbser(2023)가 점성·표면장력을 포함한 1차 쌍곡형 이상유동 모델에 쓴 엇갈린 격자(staggered grid) 반음해 스킴을 읽고, 그중 컬 구속조건(curl involution) 부분만 떼어내 재현한 기록이다. 왜 이 값이 0이 되어야 하는지, 왜 보통의 배치에서는 0이 아닌지, 그리고 무엇을 바꾸면 반올림 오차까지 0이 되는지를 숫자로 확인한다.

계면 장을 독립 변수로 승격시키는 이유#

확산계면(diffuse interface) 모델은 두 상의 경계를 색함수 $c$의 완만한 전이로 표현한다. 표면장력을 계산하려면 $\nabla c$와 그 발산이 필요하다. 문제는 미분 횟수다. 곡률은 $c$의 2차 미분이므로, 계면이 몇 격자 폭밖에 안 되면 잡음이 그대로 증폭된다.

1차 쌍곡형(first order hyperbolic) 재정식화는 여기서 미분을 한 단계 낮춘다. $\mathbf{b} = \nabla c$를 아예 독립 변수로 올리고, $\mathbf{b}$ 자신의 이송 방정식을 푼다. $c$의 이송 방정식 $\partial_t c + \mathbf{v}\cdot\nabla c = 0$의 양변에 $\nabla$를 씌우면 이렇게 된다.

$b_k$는 계면 장의 $k$ 성분, $v_j$는 속도, $v_j b_j = \mathbf{v}\cdot\mathbf{b}$는 스칼라 하나다. 이제 전체 시스템이 1차 미분만 포함한다. 고유값이 정의되고, 쌍곡형 솔버를 그대로 쓸 수 있다. 구속이 깨지면 아무 일도 없어야 할 압력도 같이 튄다. 비보존항 하나 때문에 압력이 21% 올라간 자리가 같은 종류의 장부다.

값을 하나 치른다. $\mathbf{b}$는 이제 $c$와 독립적으로 저장된다. 둘이 어긋날 수 있다.

구속조건은 초기조건이 아니라 방정식의 일부다#

$\mathbf{b}$가 어떤 스칼라의 gradient라면, 회전은 항상 0이다.

이 관계를 구속조건(involution)이라고 부른다. 이송 방정식 자체가 이 성질을 보존한다. $\partial_t (\nabla \times \mathbf{b}) = -\nabla \times \nabla(v_j b_j) = \mathbf{0}$이기 때문이다. 연속 방정식 수준에서는 $t=0$에 0이면 영원히 0이다.

그래서 흔히 초기조건 문제로 취급된다. 이것이 정확히 함정이다. 이산화된 방정식에서는 $\nabla_h \times \nabla_h$가 0이라는 보장이 없다. 보장이 없으면 매 스텝 조금씩 새고, 400스텝 뒤에는 처음 이야기한 0.544가 된다.

컬이 0이 아니면 무엇이 깨지는가. Stokes 정리가 답을 준다. 닫힌 고리를 따라 적분한 값이 그 고리가 감싼 컬의 면적분과 같다. 컬이 있으면 $\mathbf{b}$는 어떤 색함수의 gradient가 아니다. 즉 저장된 $\mathbf{b}$에 대응하는 계면이 세상에 존재하지 않는다.

이산 연산자가 Schwarz 항등식을 공유하는 자리#

논문의 처방은 단순하다. 이산 gradient와 이산 curl이 같은 네 점을 쓰게 배치하는 것이다.

색함수를 셀 중심에 두고, $b_x$는 $x$면 $(i+\frac12, j)$에, $b_y$는 $y$면 $(i, j+\frac12)$에 둔다. 컬은 노드 $(i+\frac12, j+\frac12)$에서 계산한다. 그러면 스칼라 $\psi$에 대해 이렇게 전개된다.

두 분자는 같은 네 항을 순서만 바꿔 적은 것이다. 차이는 항등적으로 0이다. 이것이 교차미분의 Schwarz 항등식(mixed partials commute)을 이산 수준에서 그대로 옮긴 형태다. 근사가 아니다. 부동소수점에서도 같은 네 수를 더하고 빼므로 반올림까지만 남는다.

이 성질이 시간 적분으로 옮겨가는 논리는 한 줄이다. 매 스텝의 갱신량이 $-\Delta t\, \nabla_h \psi$ 꼴이면, 컬을 씌운 순간 사라진다. 즉 $\nabla_h \times \mathbf{b}$는 스텝 수와 무관하게 초기값을 유지한다. 초기값이 0이면 끝까지 0이다.

Python으로 같은 전단층에 두 이산화를 걸었다#

$64^2$ 주기 격자에 동그란 계면 하나를 놓고, tanh 전단층으로 밀었다. 왼쪽은 위의 엇갈린 연산자 쌍, 오른쪽은 셀 중심에 $\mathbf{b}$를 모두 두고 성분별 Rusanov 플럭스를 쓴 배치다. 둘 다 같은 $\Delta t$, 같은 속도장이다.

import numpy as np
 
N, L = 64, 1.0
h = L / N
xc = (np.arange(N) + 0.5) * h
Xc, Yc = np.meshgrid(xc, xc, indexing="ij")
 
def shear_velocity(amp=1.0, delta=0.12):
    """Kelvin-Helmholtz 형태: tanh 전단층 + 사인 교란"""
    return amp * np.tanh((Yc - 0.5) / delta), 0.35 * amp * np.sin(2 * np.pi * Xc)
 
def interface_potential(r0=0.22, w=0.05):
    r = np.sqrt((Xc - 0.5) ** 2 + (Yc - 0.5) ** 2)
    return 0.5 * (1.0 + np.tanh((r0 - r) / w))
 
def grad_stagger(phi):
    """b_x 는 x면 (i+1/2, j), b_y 는 y면 (i, j+1/2)"""
    return (np.roll(phi, -1, 0) - phi) / h, (np.roll(phi, -1, 1) - phi) / h
 
def curl_node(bx, by):
    """노드 (i+1/2, j+1/2)의 컬 — grad_stagger 와 짝이 맞는 연산자"""
    return (np.roll(by, -1, 0) - by) / h - (np.roll(bx, -1, 1) - bx) / h
 
def curl_centered(bx, by):
    """충돌 배치의 컬: 셀 중심 성분의 중심차분"""
    return ((np.roll(by, -1, 0) - np.roll(by, 1, 0))
            - (np.roll(bx, -1, 1) - np.roll(bx, 1, 1))) / (2 * h)
 
def psi_cell(bx, by, vx, vy, upwind):
    """스칼라 하나 psi = v . b 를 셀 중심에 모은다"""
    if upwind:
        bxc = np.where(vx > 0, np.roll(bx, 1, 0), bx)
        byc = np.where(vy > 0, np.roll(by, 1, 1), by)
    else:
        bxc = 0.5 * (bx + np.roll(bx, 1, 0))
        byc = 0.5 * (by + np.roll(by, 1, 1))
    return vx * bxc + vy * byc
 
def update_compatible(bx, by, vx, vy, dt, upwind=True):
    """b <- b - dt * grad_h(psi): 더해지는 항이 전부 이산 gradient 다"""
    psi = psi_cell(bx, by, vx, vy, upwind)
    return (bx - dt * (np.roll(psi, -1, 0) - psi) / h,
            by - dt * (np.roll(psi, -1, 1) - psi) / h)
 
def update_collocated(bx, by, vx, vy, dt, lam):
    """성분별 Rusanov: 소산항은 어떤 스칼라의 gradient 도 아니다"""
    psi = vx * bx + vy * by
    fx = 0.5 * (psi + np.roll(psi, -1, 0)) - 0.5 * lam * (np.roll(bx, -1, 0) - bx)
    fy = 0.5 * (psi + np.roll(psi, -1, 1)) - 0.5 * lam * (np.roll(by, -1, 1) - by)
    return (bx - dt * (fx - np.roll(fx, 1, 0)) / h,
            by - dt * (fy - np.roll(fy, 1, 1)) / h)
 
def path_gap(bx, by):
    """같은 셀까지 두 경로로 b 를 적분한다. 차이가 감싼 컬이다"""
    i1 = j1 = N // 2
    pa = bx[:i1, 0].sum() * h + by[i1, :j1].sum() * h   # x 먼저, 그다음 y
    pb = by[0, :j1].sum() * h + bx[:i1, j1].sum() * h   # y 먼저, 그다음 x
    return pa - pb
 
vx, vy = shear_velocity()
phi0 = interface_potential()
lam = float(np.max(np.abs(vx) + np.abs(vy)))
dt = 0.4 * h / lam
bs = list(grad_stagger(phi0))
bc = [np.gradient(phi0, h, axis=0, edge_order=2),
      np.gradient(phi0, h, axis=1, edge_order=2)]
 
print(f"N={N}  h={h:.4e}  dt={dt:.4e}  lambda={lam:.3f}")
print(f"{'t':>6}{'|b|max stag':>13}{'|curl|max stag':>16}"
      f"{'|b|max coll':>13}{'|curl|max coll':>16}")
for k in range(5):
    if k:
        for _ in range(100):
            bs = update_compatible(*bs, vx, vy, dt)
            bc = update_collocated(*bc, vx, vy, dt, lam)
    print(f"{k*100*dt:>6.3f}{np.max(np.hypot(*bs)):>13.3f}"
          f"{np.max(np.abs(curl_node(*bs))):>16.3e}"
          f"{np.max(np.hypot(*bc)):>13.3f}"
          f"{np.max(np.abs(curl_centered(*bc))):>16.3e}")
 
print(f"\npath-integral gap of phi at (0.5, 0.5):"
      f"  staggered {path_gap(*bs):+.3e}   collocated {path_gap(*bc):+.3e}")
 
bn = list(grad_stagger(phi0))
for _ in range(400):
    bn = update_compatible(*bn, vx, vy, dt, upwind=False)
print(f"staggered with NO upwinding: |b|max={np.max(np.hypot(*bn)):.3e} (unstable)"
      f"  |curl|max={np.max(np.abs(curl_node(*bn))):.3e}")

N=64  h=1.5625e-02  dt=4.6330e-03  lambda=1.349
     t  |b|max stag  |curl|max stag  |b|max coll  |curl|max coll
 0.000       10.167       0.000e+00        9.818       3.286e-04
 0.463       12.318       1.421e-12        6.538       5.545e+01
 0.927       18.063       1.748e-12        4.703       1.687e+02
 1.390        9.854       1.705e-12        4.124       2.269e+02
 1.853        5.746       2.132e-12        3.342       1.998e+02
 
path-integral gap of phi at (0.5, 0.5):  staggered -3.331e-16   collocated -5.442e-01
staggered with NO upwinding: |b|max=3.622e+03 (unstable)  |curl|max=1.091e-10

$\kappa$ 슬라이더를 올리면 두 경로의 막대가 갈라지고, 색칠된 면의 컬 적분이 그 차와 같이 커진다.

0으로 내리면 두 막대가 같은 숫자에 멈춘다. 조금만 올려도 경로 차가 그 면적분과 맞는다.

엇갈린 쪽의 컬은 $10^{-12}$에 머문다. $1/h^2 \approx 4096$이 곱해진 배정밀도 반올림이 그 크기다. 충돌 쪽은 0.46 시간만에 55, 1.39에서 227까지 올라간다. $|\mathbf{b}|$의 최대값이 4 남짓이니 컬은 장 자체보다 50배 이상 크다. 남아 있는 것은 계면이 아니라 회전 잡음이다.

$t=1.390$에서 충돌 배치의 컬 최댓값은 226.9다. 같은 시각 $|\mathbf{b}|$ 최댓값은 4.124다. 둘을 나누면 50을 넘긴다. 장보다 회전이 크다. 이 컬을 $(0.5, 0.5)$까지 두 경로로 적분한 값이 경로 차 $-5.442\times 10^{-1}$이다. 색함수는 0과 1 사이인데, 경로 차가 그 폭의 절반을 넘는다.

시간 칸도 표에 그대로 있다. $dt=4.6330\times 10^{-3}$이다. 100스텝이 0.463이고 400스텝이 1.853이다. 출력의 시각 열과 같다. 그 사이 충돌 컬은 $3.286\times 10^{-4}$에서 55.45, 168.7, 226.9로 올랐다가 1.853에서 199.8이 된다. $|\mathbf{b}|$가 3.342까지 깎여 컬의 절대값도 조금 줄었다. 컬이 사라진 것은 아니다. 엇갈린 쪽은 같은 400스텝 동안 $1.421\times 10^{-12}$에서 $2.132\times 10^{-12}$다. $1/h^2 \approx 4096$이 배정도 반올림에 곱해진 크기와 맞다.

상류를 끈 줄은 다른 실패다. $|\mathbf{b}|$는 $3.622\times 10^{3}$까지 커진다. 초기 10.167의 약 360배다. 컬은 $1.091\times 10^{-10}$에 남는다. 표의 $\lambda=1.349$와 $h=1.5625\times 10^{-2}$를 곱한 크기가, 갱신이 스칼라 하나의 기울기가 아닐 때 한 스텝에 새는 비율이다. 400스텝이면 그 누적이 경로 차 0.544까지 간다. 항등식은 크기를 잡아 주지 않는다. 소산은 따로 넣어야 한다.

소산을 키워도 컬이 늘지 않는 조건#

두 번째 줄의 update_compatible은 상류차분(upwind)을 쓴다. 소산이 들어간다는 뜻이다. 그런데 컬은 여전히 $10^{-12}$다. 여기가 이 스킴의 실질적인 이점이다.

이유는 상류차분이 들어간 자리에 있다. 상류 판정은 스칼라 $\psi = \mathbf{v}\cdot\mathbf{b}$를 셀 중심에 모을 때만 쓰인다. 그다음 갱신은 여전히 $-\Delta t\, \nabla_h \psi$ 한 덩어리다. $\psi$를 어떻게 만들든 gradient의 컬은 0이다.

충돌 배치는 반대다. Rusanov 소산항 $-\frac{\lambda}{2}(b_{x,i+1} - b_{x,i})$는 $b_x$에는 $x$방향으로, $b_y$에는 $y$방향으로만 붙는다. 두 성분이 서로 다른 연산자를 받는다. 이 쌍은 어떤 스칼라의 gradient가 아니므로, 컬을 씌우면 $\frac{\lambda h}{2}(\partial_x \partial_y^2 b_y - \partial_y \partial_x^2 b_x)$가 남는다.

제약은 정확도가 아니라 형태다. 갱신량 전체가 스칼라 하나의 이산 gradient이기만 하면 리미터든 고차 재구성이든 얹어도 된다. 격자가 표현할 수 있는 항을 먼저 세고, 스킴은 그다음에 고른다.

컬이 0이어도 계산은 터진다#

마지막 출력 줄이 이 이야기를 한다. 상류차분을 끄고 중심평균만 쓰면 $|\mathbf{b}|$가 $3.6\times 10^3$까지 자란다. 초기값이 10이었으니 360배다. 계산은 발산했다. 그런데 컬은 $1.09\times 10^{-10}$다. 여전히 반올림 수준이다.

구속조건 보존과 안정성은 별개의 성질이다. 엇갈린 배치는 "이 장이 여전히 어떤 계면의 gradient다"를 보장할 뿐, 그 계면이 유한한 크기에 머문다고 약속하지 않는다. 소산은 따로 넣어야 한다. 둘을 섞어 생각하면, 발산한 계산을 보고 격자 배치를 의심하는 잘못된 진단에 이른다.

표면장력에서 이 성질이 값을 바꾸는 자리#

논문이 다루는 시스템은 $\mathbf{b}$로부터 표면장력 응력을 직접 만든다. 응력 텐서는 대략 $\sigma(|\mathbf{b}|\,\mathbf{I} - \mathbf{b}\otimes\mathbf{b}/|\mathbf{b}|)$ 꼴이다. $\mathbf{b}$에 회전 성분이 섞이면 이 응력이 존재하지 않는 계면의 장력을 만든다.

증상은 정적 기포에서 먼저 보인다. 이론상 속도가 0이어야 하는 문제에서 기생류(spurious current)가 자란다. 컬 성분은 곡률 계산을 통해 응력에 들어가고, 응력은 속도를 만들고, 속도는 다시 $\mathbf{b}$를 회전시킨다. 되먹임 고리가 닫히면 진폭이 지수적으로 커진다.

논문은 같은 격자 위에서 압력을 반음해로 푼다. 음향파만 음해로 넘기고 나머지는 명시적으로 둔다. 엇갈린 배치가 여기서 두 번 일한다. 컬 항등식을 주고, 압력 방정식을 대칭양정치로 만들어 켤레구배법(CG)으로 풀 수 있게 한다. 배치 하나가 두 성질을 동시에 사준다.

컬이 0인 채로 터지면 격자를 의심하지 마라#

$b_x$와 $b_y$가 다른 면에 살아서 $\mathbf{v}\cdot\mathbf{b}$를 만들 때마다 평균이 들어간다. 코드의 psi_cell이 그 평균이다. 평균은 고주파를 깎는다. 계면이 3~4칸이면 이 감쇠를 빼 놓고 배치를 고르면 계면 폭이 먼저 죽는다. 면과 노드의 경계 스텐실도 따로다. 비정렬 격자로 옮기면 호환 연산자 쌍을 다시 짜야 해서, 논문이 직교 격자에 머문다.

판단은 둘로 갈린다. 계면 장이나 변형장을 독립 변수로 두면 컬 항등식은 선택이 아니다. 매 스텝 새고, 새는 속도는 $\lambda h$에 비례한다. $\mathbf{b}$를 매번 $c$에서 다시 구하는 모델이면 이 누수는 생기지 않는다. 어느 쪽이든 path_gap 열 줄로, 저장한 벡터장에 대응하는 스칼라가 있는지 한 번은 센다. 그 격자에 압력을 음해로 얹는 비용은 시간스텝을 1001배 키운 반음해 압력과 한 장부다.