Con el área curva exacta de cada celda, ∂x/∂x dejó de valer 1 — por qué el área por producto cruz de diagonales encaja con Green–Gauss

Un área más exacta dio una respuesta peor#

El gradiente de $\phi = x$ vale 1. Una malla curva no cambia eso. Al integrar el área de celda sobre las curvas reales en un FVM curvilíneo, el gradiente se desvió un 0.19%. El jacobiano en el centro de la celda abre la brecha hasta 0.39%. El producto cruz de las dos diagonales, con cuatro vértices, se queda en $10^{-14}$.

El punto de partida es un problema de tarea de CFD de posgrado. Hay que hallar el $\Delta V$ 2D con el producto cruz de dos diagonales. Cinco líneas a mano. La tarea no pregunta por qué esas líneas superan a la integral curva. La razón es una identidad del teorema de la divergencia.

Al arrastrar los vértices A–D, o al pulsar sweep C, el área con signo de la esquina superior izquierda se mueve junto con el paralelogramo morado.

Cuando C pasa a la izquierda de D, los colores de los triángulos se separan y el signo del área cae a negativo. Ese signo es la misma clase de señal que $-0.500$ en $x=-2$, al final de la salida.

La respuesta de la tarea: dos diagonales dan el área del cuadrilátero#

La celda $(i,j)$ del espacio computacional se convierte en un cuadrilátero ABCD en el espacio físico. Sus vértices están en índices semienteros $(i\pm\tfrac12, j\pm\tfrac12)$. En 2D se toma profundidad unitaria, así que el volumen coincide con el área.

Las diagonales $\overline{BD}$ y $\overline{CA}$ se tratan como vectores. Al trasladarlas a un origen común forman el paralelogramo EFGH. Su área es la magnitud del producto cruz $|\mathbf{BD}\times\mathbf{CA}|$. Ese paralelogramo mide exactamente el doble que el cuadrilátero original.

El primer par de paréntesis es la componente $x$ de una diagonal por la componente $y$ de la otra. El segundo par es la combinación opuesta. La diferencia es el producto cruz 2D (el área con signo del paralelogramo que generan dos vectores).

Por qué el producto cruz de diagonales no es una aproximación#

Al cortar el cuadrilátero por la diagonal AC se obtienen los triángulos ABC y ACD. El área de cada uno es la mitad del producto cruz de sus vectores de lado. Al sumarlas, se agrupan así:

Como $\mathbf{AD}-\mathbf{AB}=\mathbf{BD}$, se llega a la última forma. La identidad vale para cuadriláteros convexos y cóncavos. Por eso la fórmula de diagonales es exactamente el área del polígono que une los vértices con rectas. Coincide con la fórmula del cordón de zapato (shoelace, la suma cruzada de coordenadas de vértices para el área de un polígono).

Pero esa es el área de un polígono de lados rectos. Si las líneas de malla son realmente curvas, el área real es $\iint J\,d\xi\,d\eta$. Aquí $J$ es el jacobiano (el factor de estiramiento de área entre coordenadas computacionales y físicas). Ambas difieren en $O(h^2)$. ¿Es entonces mejor el área verdadera? Todo lo contrario.

Lo que el gradiente de Green–Gauss exige al volumen#

El gradiente de Green–Gauss (que usa el teorema de la divergencia para convertir el gradiente medio de celda en una suma sobre caras) se escribe así:

$\phi_f$ es el valor en el punto medio de la cara, $\mathbf{S}_f$ es la normal de cara por su longitud y $V_P$ es el volumen de celda. En un código FVM las caras son segmentos rectos entre vértices. Ninguna cara sigue la línea de malla curva.

Con $\phi=x$: sobre un lado recto, $x$ es lineal, así que su valor en el punto medio es exactamente el promedio del lado. El numerador es entonces el teorema de la divergencia aplicado al polígono de lados rectos.

$\Omega_\text{poly}$ es el polígono que une los vértices con rectas. Así, el gradiente vale $|\Omega_\text{poly}|/V_P$. Para que sea 1, $V_P$ debe ser el área del polígono, es decir, el área por producto cruz de diagonales. Con el área curva verdadera, numerador y denominador miden figuras distintas.

El mismo desajuste aparece al discretizar los términos métricos de transformación de coordenadas curvilíneas y métricas de malla. Los vectores de cara y los volúmenes deben salir de los mismos datos de vértices.

Código: mismas caras, tres volúmenes#

La malla ondulada es $x=\xi+A\sin 2\pi\xi\,\sin 2\pi\eta$, $y=\eta+A\sin 2\pi\xi\,\sin 2\pi\eta$ con $A=0.06$. La frontera queda recta y solo se curvan las líneas interiores. Los flujos de cara son idénticos en los tres casos; solo cambia el volumen: (1) el producto cruz de diagonales, (2) el $J\cdot h^2$ analítico en el centro de celda, (3) el área curva verdadera con cuadratura de Gauss 4×4.

import math
 
A = 0.06  # amplitud de la ondulación de la malla
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): producto cruz de las dos diagonales
    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 analítico en el centro de celda por 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):
    # área de la celda curva: Gauss-Legendre 4x4 sobre J
    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):
    # d(phi)/dx para phi = x; caras rectas, antihorario
    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]   # normal exterior (x) por longitud de cara
        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

Al cambiar el volumen entre ½ (BD × CA) y el $J$ en el centro de celda, el flujo de cara se queda fijo y solo se parten los colores.

Las curvas blancas son las líneas de malla del mapeo. Los cuadriláteros de color son las celdas FV de lados rectos. Con ½ (BD × CA) todas las celdas quedan verdes. Con $J$, las celdas más curvadas se ponen rojas. Subir N apaga el rojo. No lo quita.

Entre 1.1e-14 y 3.9e-3#

Con $N=32$, el área por diagonales da un error de $1.1\times10^{-14}$, nivel de redondeo. El $J$ en el centro da $3.881\times10^{-3}$ y el área curva verdadera, $1.935\times10^{-3}$. Ambos caen a la cuarta parte al refinar la malla al doble. Es orden 2.00.

El área total vale 1 con los tres métodos. La suma global es correcta. El reparto por celda no lo es. Una verificación de conservación global no detecta este error.

El 0.19% en $N=32$ no es otra medición. Es el error del área curva $1.935\times10^{-3}$, leído como porcentaje. En esa misma fila, el $J$ de centro $3.881\times10^{-3}$ es el 0.39% del inicio. La casilla de las diagonales es $1.1\times10^{-14}$. Las caras coinciden y el numerador también. Solo cambia el denominador $V_P$.

La columna de $J$, leída hacia abajo, empieza en $6.032\times10^{-2}$ para $N=8$. Siguen $1.543\times10^{-2}$, $3.881\times10^{-3}$, $9.716\times10^{-4}$, $2.430\times10^{-4}$ en $N=16,32,64,128$. Los órdenes impresos son 1.97, 1.99, 2.00, 2.00. De $N=8$ a $N=128$ la malla se afina 16 veces y el error baja unas 248 veces. El segundo orden marca $4^{4}=256$. La distancia entre 248 y 256 cabe en ese primer 1.97. Después de $N=8$, cada paso es casi exactamente un factor de cuatro.

La columna del área curva baja en paralelo: $2.896\times10^{-2}$, $7.638\times10^{-3}$, $1.935\times10^{-3}$, $4.855\times10^{-4}$, $1.215\times10^{-4}$. Cada fila queda cerca de la mitad del error del $J$ central. En 8 celdas, $2.896/6.032\approx0.48$. En 128, $1.215/2.430=0.50$. Una cuadratura de Gauss 4×4 supera una muestra en el centro. El numerador sigue midiendo el polígono recto $|\Omega_\text{poly}|$, así que el gradiente $|\Omega_\text{poly}|/V_P$ se aparta de 1. Ambos totales imprimen 1.000000000000000. El exceso y el déficit por celda se cancelan. El $\partial x/\partial x$ de una celda no ve esa cancelación.

Un error $O(h^2)$ se deja pasar con facilidad dentro de un esquema de segundo orden. Si un campo lineal no se reproduce con exactitud, la premisa de la reconstrucción de segundo orden se tambalea. Los limitadores también reaccionan a este gradiente falso. Se suma a la oblicuidad vista en gradientes de Green–Gauss frente a mínimos cuadrados.

La señal que da esta fórmula en una celda torcida#

Las últimas cuatro líneas de la salida salen de arrastrar el vértice C hacia la izquierda. Con C en $x=+1$ el área es 1.000, y en $x=0$ es 0.500. En $x=-1$ es 0.000, y en $x=-2$ es $-0.500$.

Cuando C pasa a D, los lados BC y DA se cruzan. La celda se tuerce en forma de pajarita (bow-tie). La fórmula de diagonales devuelve entonces la suma con signo de las áreas de los dos triángulos. Un triángulo se vuelve negativo, así que la suma cae a cero o por debajo.

No es un defecto. Es una señal. El signo de $\Delta V$, revisado en todas las celdas justo después de generar la malla, marca las torcidas. El valor absoluto apaga la señal. La celda torcida entra al cálculo como una celda normal de volumen positivo.

Si ∂x/∂x no vale 1, ese volumen no entra en la actualización#

El $J$ analítico no sobra. Sigue sirviendo como indicador de calidad de malla y como estimador de error. El $V_P$ del denominador de la actualización FVM tiene que salir de los mismos vértices que los vectores de cara $\mathbf{S}_f$. Meter la integral curva $\iint J\,d\xi\,d\eta$ tal cual hace que el denominador mire otra figura.

La prueba es corta. Se colocan $\phi=x$ y $\phi=y$ en una malla curvada y se imprime el gradiente de Green–Gauss. Si $\partial x/\partial x$ no vale 1 a nivel de redondeo, ese volumen no entra en la actualización. Las cinco líneas de diagonales de la tarea son el camino corto para pasar la prueba. El mismo tipo de desajuste, un área que rompe el balance y deja mover un estado en reposo, aparece en el FVM polar well-balanced.