곡선 셀의 참 면적을 넣었더니 ∂x/∂x가 1에서 어긋났다 — 대각선 외적 면적이 Green–Gauss와 맞물리는 이유

더 정확한 면적이 더 나쁜 답을 냈다#

$\phi = x$의 기울기는 1이다. 격자가 휘어도 그렇다. 곡선격자 FVM에 셀 면적을 곡선 그대로 적분해 넣었더니 0.19%가 빠졌다. 셀 중심 야코비안은 0.39%로 더 벌어졌다. 꼭짓점 네 개의 대각선 외적은 $10^{-14}$에서 멈춘다.

출발점은 대학원 CFD 숙제다. 2차원 $\Delta V$를 대각선 두 개의 외적으로 구하라는 것. 손으로 다섯 줄. 숙제는 그 다섯 줄이 왜 곡선 적분보다 나은지 묻지 않는다. 이유는 발산 정리 하나다.

꼭짓점 A–D를 끌거나 sweep C를 누르면, 왼쪽 위 부호 있는 면적과 보라색 평행사변형이 같이 움직인다.

C가 D 왼쪽으로 넘어가면 삼각형 색이 갈라지고 면적 부호가 음수로 떨어진다. 그 부호가 출력 맨 아래, $x=-2$의 $-0.500$과 같은 종류의 신호다.

숙제의 답: 사각형 면적은 대각선 두 개로 끝난다#

계산 공간의 셀 $(i,j)$는 실제 공간에서 꼭짓점 네 개를 가진 사각형 ABCD가 된다. 꼭짓점은 반정수 인덱스 $(i\pm\tfrac12, j\pm\tfrac12)$에 있다. 2차원이므로 높이를 1로 두면 부피는 곧 면적이다.

대각선 $\overline{BD}$와 $\overline{CA}$를 벡터로 보자. 두 벡터를 평행이동해 시작점을 맞추면 평행사변형 EFGH가 생긴다. 그 넓이는 외적의 크기 $|\mathbf{BD}\times\mathbf{CA}|$이다. 그리고 이 평행사변형은 원래 사각형 넓이의 정확히 두 배다.

첫 괄호 쌍은 한 대각선의 $x$ 성분과 다른 대각선의 $y$ 성분이다. 둘째 쌍은 그 반대 조합이다. 차이가 2차원 외적(cross product, 두 벡터가 만드는 평행사변형의 부호 있는 넓이)이다.

왜 대각선 외적은 근사가 아닌가#

사각형을 대각선 AC로 자르면 삼각형 ABC와 ACD가 된다. 각 삼각형 넓이는 변 벡터 외적의 절반이다. 둘을 더하면 이렇게 묶인다.

$\mathbf{AD}-\mathbf{AB}=\mathbf{BD}$이므로 마지막 식이 나온다. 볼록이든 오목이든 이 등식은 그대로다. 그러니 대각선 공식은 꼭짓점을 직선으로 이은 다각형의 넓이와 정확히 같다. 신발끈 공식(shoelace, 꼭짓점 좌표를 엇갈려 곱해 다각형 넓이를 구하는 공식)과도 같은 값이다.

단, 이것은 "직선 변" 다각형의 넓이다. 격자선이 실제로 곡선이면 참 면적은 $\iint J\,d\xi\,d\eta$이다. $J$는 야코비안(계산 좌표 → 실제 좌표 변환의 면적 확대율)이다. 두 값은 $O(h^2)$만큼 다르다. 그럼 참 면적이 더 좋은 선택일까. 답은 반대다.

Green–Gauss 기울기가 부피에 요구하는 것#

Green–Gauss 기울기(발산 정리로 셀 평균 기울기를 면 적분 합으로 바꾸는 방법)는 이렇게 쓴다.

$\phi_f$는 면 중점 값, $\mathbf{S}_f$는 면 법선 벡터에 면 길이를 곱한 것, $V_P$는 셀 부피다. FVM 코드에서 면은 꼭짓점을 잇는 직선이다. 곡선 격자선을 따라가는 면은 없다.

$\phi=x$를 넣어 보자. 직선 변 위에서 $x$는 선형이므로 중점 값이 변 평균과 정확히 같다. 그러면 분자는 직선 다각형에 대한 발산 정리 그 자체다.

$\Omega_\text{poly}$는 꼭짓점을 직선으로 이은 다각형이다. 결국 기울기는 $|\Omega_\text{poly}|/V_P$가 된다. 이 값이 1이려면 $V_P$는 다각형 넓이여야 한다. 곧 대각선 외적 면적이다. 참 곡선 면적을 넣으면 분자와 분모가 서로 다른 도형을 재게 된다.

이 불일치는 곡선 좌표 변환과 격자 메트릭에서 다룬 메트릭 항을 이산화할 때도 같은 형태로 나타난다. 면 벡터와 부피를 같은 꼭짓점 데이터에서 만들어야 한다는 조건이다.

코드: 같은 면, 세 가지 부피#

물결 격자 $x=\xi+A\sin 2\pi\xi\,\sin 2\pi\eta$, $y=\eta+A\sin 2\pi\xi\,\sin 2\pi\eta$를 쓴다. $A=0.06$이다. 경계는 직선으로 남고 내부 격자선만 휜다. 면 플럭스는 셋 다 같고 부피만 바꾼다. (1) 대각선 외적, (2) 셀 중심 해석 $J\cdot h^2$, (3) 4×4 가우스 적분으로 구한 참 곡선 면적이다.

import math
 
A = 0.06  # 격자 물결 진폭
TP = 2 * math.pi
 
def mapping(xi, eta):
    w = A * math.sin(TP * xi) * math.sin(TP * eta)
    return xi + w, eta + w
 
def jac(xi, eta):
    sx, cx = math.sin(TP * xi), math.cos(TP * xi)
    se, ce = math.sin(TP * eta), math.cos(TP * eta)
    wx, we = A * TP * cx * se, A * TP * sx * ce
    return (1 + wx) * (1 + we) - we * wx
 
def wavy_nodes(n):
    s = [i / n for i in range(n + 1)]
    x = [[mapping(si, sj)[0] for sj in s] for si in s]
    y = [[mapping(si, sj)[1] for sj in s] for si in s]
    return s, x, y
 
def diag_area(x, y, i, j):
    # 1/2 (BD x CA): 두 대각선의 외적
    d1x = x[i + 1][j + 1] - x[i][j]; d1y = y[i + 1][j + 1] - y[i][j]
    d2x = x[i][j + 1] - x[i + 1][j]; d2y = y[i][j + 1] - y[i + 1][j]
    return 0.5 * (d1x * d2y - d1y * d2x)
 
def jacobian_area(s, i, j):
    # 셀 중심의 해석적 J에 h^2를 곱함
    h = s[1] - s[0]
    return jac(s[i] + h / 2, s[j] + h / 2) * h * h
 
GL = [(-0.8611363115940526, 0.3478548451374538), (-0.3399810435848563, 0.6521451548625461),
      (0.3399810435848563, 0.6521451548625461), (0.8611363115940526, 0.3478548451374538)]
 
def exact_cell_area(s, i, j):
    # 곡선 셀의 면적: J에 4x4 가우스-르장드르 적분
    h = s[1] - s[0]; area = 0.0
    for a, wa in GL:
        for b, wb in GL:
            area += wa * wb * jac(s[i] + h * (1 + a) / 2, s[j] + h * (1 + b) / 2)
    return area * h * h / 4
 
def green_gauss_dx(x, y, i, j, vol):
    # phi = x의 d(phi)/dx; 직선 면, 반시계 방향
    p = [(i, j), (i + 1, j), (i + 1, j + 1), (i, j + 1)]
    flux = 0.0
    for k in range(4):
        (ia, ja), (ib, jb) = p[k], p[(k + 1) % 4]
        nx = y[ib][jb] - y[ia][ja]   # 바깥 법선(x 성분) x 면 길이
        flux += 0.5 * (x[ia][ja] + x[ib][jb]) * nx
    return flux / vol
 
print("  N   err(diag)  err(J centre)  err(exact curved)  order(J)")
prev = None
for n in (8, 16, 32, 64, 128):
    s, x, y = wavy_nodes(n)
    e = [0.0, 0.0, 0.0]
    for i in range(n):
        for j in range(n):
            vols = (diag_area(x, y, i, j), jacobian_area(s, i, j), exact_cell_area(s, i, j))
            for k, v in enumerate(vols):
                e[k] = max(e[k], abs(green_gauss_dx(x, y, i, j, v) - 1.0))
    rate = "   -" if prev is None else f"  {math.log2(prev / e[1]):.2f}"
    print(f"{n:3d}   {e[0]:.1e}    {e[1]:.3e}      {e[2]:.3e}     {rate}")
    prev = e[1]
 
s, x, y = wavy_nodes(32)
tot_d = sum(diag_area(x, y, i, j) for i in range(32) for j in range(32))
tot_e = sum(exact_cell_area(s, i, j) for i in range(32) for j in range(32))
print(f"total area  diag {tot_d:.15f}  exact {tot_e:.15f}")
 
def bowtie_area(pts):
    (ax, ay), (bx, by), (cx, cy), (dx, dy) = pts
    return 0.5 * ((cx - ax) * (dy - by) - (cy - ay) * (dx - bx))
 
for t in (0.0, 1.0, 2.0, 3.0):
    a = bowtie_area([(0, 0), (1, 0), (1 - t, 1), (0, 1)])
    print(f"C at x={1 - t:+.0f}: signed area {a:+.3f}")

  N   err(diag)  err(J centre)  err(exact curved)  order(J)
  8   2.3e-15    6.032e-02      2.896e-02        -
 16   3.6e-15    1.543e-02      7.638e-03       1.97
 32   1.1e-14    3.881e-03      1.935e-03       1.99
 64   3.4e-14    9.716e-04      4.855e-04       2.00
128   6.8e-14    2.430e-04      1.215e-04       2.00
total area  diag 1.000000000000000  exact 1.000000000000000
C at x=+1: signed area +1.000
C at x=+0: signed area +0.500
C at x=-1: signed area +0.000
C at x=-2: signed area -0.500

부피 방법을 ½ (BD × CA)와 셀 중심 $J$ 사이에서 바꾸면, 같은 면 플럭스인데 셀 색만 갈린다.

흰 곡선은 매핑의 격자선이고, 색 사각형은 직선 변 FV 셀이다. ½ (BD × CA)를 고르면 전 셀이 녹색이다. $J$ 기반은 많이 휜 셀이 붉다. N을 올리면 붉은색은 옅어질 뿐 사라지지 않는다.

1.1e-14와 3.9e-3 사이#

$N=32$에서 대각선 면적은 오차 $1.1\times10^{-14}$다. 반올림 수준이다. 셀 중심 $J$는 $3.881\times10^{-3}$, 참 곡선 면적은 $1.935\times10^{-3}$이다. 둘 다 격자를 두 배 조이면 4분의 1이 된다. 차수 2.00이다.

전체 면적 합은 세 방법 모두 1이다. 총량은 맞고 셀별 분배가 틀리다. 전역 보존 검사로는 이 오류가 안 보인다.

$N=32$의 0.19%는 따로 잰 값이 아니다. 참 곡선 면적 오차 $1.935\times10^{-3}$을 퍼센트로 읽은 것이다. 같은 줄의 셀 중심 $J$인 $3.881\times10^{-3}$이 도입의 0.39%다. 대각선 칸은 $1.1\times10^{-14}$다. 면은 같고 분자도 같다. 분모 $V_P$만 셋이다.

$J$ 열을 아래로 읽는다. $N=8$은 $6.032\times10^{-2}$다. 이어서 $N=16,32,64,128$에서 $1.543\times10^{-2}$, $3.881\times10^{-3}$, $9.716\times10^{-4}$, $2.430\times10^{-4}$다. 인쇄된 차수는 1.97, 1.99, 2.00, 2.00이다. $N$이 8에서 128로 16배가 되는 동안 오차는 약 248배가 줄었다. 2차의 기준 $4^{4}=256$과 248의 차이는 첫 칸의 1.97에 들어 있다. $N=8$ 다음부터는 거의 정확히 4분의 1이다.

참 곡선 면적 열도 같이 내려간다. $2.896\times10^{-2}$, $7.638\times10^{-3}$, $1.935\times10^{-3}$, $4.855\times10^{-4}$, $1.215\times10^{-4}$. 어느 행에서나 셀 중심 $J$의 절반 근처다. 8칸은 $2.896/6.032\approx0.48$, 128칸은 $1.215/2.430=0.50$이다. 가우스 4×4가 중심 한 점보다 나은 것은 맞다. 그래도 분자가 재는 도형은 직선 다각형 $|\Omega_\text{poly}|$라서, 기울기는 $|\Omega_\text{poly}|/V_P$에서 1을 벗어난다. 전역 합은 양쪽 다 1.000000000000000이다. 셀마다 과대와 과소가 서로 지운 결과다. 한 셀의 $\partial x/\partial x$는 그 상쇄를 보지 못한다.

$O(h^2)$라서 2차 스킴이면 괜찮다고 넘기기 쉽다. 선형장조차 정확히 재현하지 못하면 2차 재구성의 전제가 흔들린다. 리미터도 이 가짜 기울기에 반응한다. 비뚤어진 셀에서 Green–Gauss와 최소제곱이 갈라지는 비교 위에, 부피 오차가 한 겹 더 얹힌다.

비틀린 셀에서 이 공식이 내는 신호#

출력의 마지막 네 줄은 꼭짓점 C를 왼쪽으로 끌어 간 결과다. C가 $x=+1$일 때 면적은 1.000, $x=0$에서 0.500이다. $x=-1$에서는 0.000, $x=-2$에서는 $-0.500$이다.

C가 D를 지나면 변 BC와 DA가 교차한다. 셀이 나비 모양(bow-tie)으로 비틀린다. 이때 대각선 공식은 두 삼각형 넓이의 부호 있는 합을 낸다. 한쪽 삼각형이 음수가 되므로 합이 0이나 음수로 떨어진다.

이것은 결함이 아니라 신호다. 격자 생성 직후 모든 셀의 $\Delta V$ 부호를 보면 비틀린 셀이 바로 나온다. 절댓값을 씌우면 신호가 죽는다. 비틀린 셀이 양의 부피를 가진 정상 셀로 계산에 들어간다.

∂x/∂x가 1이 아니면 그 부피는 쓰지 않는다#

해석적 $J$를 버리라는 뜻은 아니다. 격자 품질이나 오차 추정에는 아직 쓸 만하다. FVM 업데이트 분모의 $V_P$만은 면 벡터 $\mathbf{S}_f$와 같은 꼭짓점에서 나와야 한다. 곡선 적분 $\iint J\,d\xi\,d\eta$를 그대로 넣으면 분모가 다른 도형을 본다.

판정은 짧다. 휜 격자에 $\phi=x$와 $\phi=y$를 넣고 Green–Gauss를 찍는다. $\partial x/\partial x$가 반올림 수준으로 1이 아니면 그 부피는 업데이트에 쓰지 않는다. 숙제의 대각선 다섯 줄이 그 검사를 통과하는 짧은 길이다. 면적이 균형식과 어긋나 정지 해가 움직이는 경우는 극좌표 FVM의 well-balanced에서도 같은 종류다.