Cuando el gradiente se tuerce en mallas deformadas — Green–Gauss y reconstrucción por mínimos cuadrados ponderados

Dos métodos para calcular el gradiente en celdas de malla no estructurada y cuál resiste mejor la distorsión

El residual empezó a aplanarse de forma extraña y luego divergió. Al revisar los logs, en una celda cerca de la pared el gradiente de temperatura apuntaba casi perpendicular al valor real. La malla mostraba que solo esa celda era anormalmente achatada. Los valores estaban bien. El método que construye el gradiente tropezó con la forma de la malla.

En el método de volúmenes finitos (FVM), el gradiente de celda alimenta el flujo difusivo, la reconstrucción de segundo orden y el limitador. El valor en el centro de la celda no basta. Hay que reunir información de los vecinos para estimar $\nabla\phi_C$ . Esta entrada compara dos métodos para hacer esa estimación — Green–Gauss y mínimos cuadrados. Se ve, con código, por qué uno se derrumba en mallas deformadas y qué corrige la ponderación por distancia.

Green–Gauss: rodear con caras y promediar#

Se aplica el teorema de la divergencia directamente a la celda. La integral de volumen se transforma en una integral de superficie.

\nabla \phi_C = \frac{1}{V_C} \sum_{f} \phi_f \, \mathbf{n}_f \, A_f

$V_C$ es el volumen de la celda, $\mathbf{n}_f$ la normal exterior a la cara, $A_f$ el área de la cara y $\phi_f$ el valor sobre la cara.

El valor de cara $\phi_f$ no se conoce directamente. Se interpola desde los valores en los centros de celda. La interpolación lineal más común es esta.

\phi_f = g_f\, \phi_C + (1 - g_f)\, \phi_N

$g_f$ es un peso definido por la razón de distancias desde la cara hasta el centro del vecino $N$ . El cálculo es barato e intuitivo. En mallas estructuradas alcanza segundo orden de exactitud.

La malla deformada sacude a Green–Gauss#

El problema se esconde en $g_f$ . La interpolación lineal supone que el centro de la cara y los dos centros de celda están sobre una misma recta. En mallas no estructuradas esa suposición casi siempre se rompe.

Si el centro de la cara se aparta de la línea que une los dos centros de celda ( $\text{skewness}$ ), el $\phi_f$ interpolado queda sesgado. Si la celda es achatada o los vecinos difieren mucho en tamaño, el sesgo crece más. Este error se acumula con las iteraciones. Al final el gradiente apunta en una dirección equivocada.

Añadir un término de corrección lo atenúa. Pero la corrección vuelve a necesitar el gradiente. Se genera un círculo. La corrección iterativa es cara y, si la distorsión es fuerte, la convergencia se vuelve lenta.

Mínimos cuadrados: ajustar un plano a los vecinos#

Se cambia el enfoque. No se pasa por las caras. Se ajusta un solo plano directamente a la distribución de valores alrededor de la celda $C$ .

Si el vecino $k$ está en el desplazamiento $\mathbf{d}_k = \mathbf{x}_k - \mathbf{x}_C$ , la aproximación lineal es $\phi_k \approx \phi_C + \nabla\phi_C \cdot \mathbf{d}_k$ . Se minimiza la suma de los cuadrados de este error sobre todos los vecinos.

\min_{\nabla\phi_C} \sum_k w_k \left( \phi_k - \phi_C - \nabla\phi_C \cdot \mathbf{d}_k \right)^2

$w_k$ es el peso por vecino. Al derivar e igualar a cero surge la ecuación normal.

\left( \sum_k w_k\, \mathbf{d}_k \mathbf{d}_k^{\mathsf T} \right) \nabla\phi_C = \sum_k w_k \left( \phi_k - \phi_C \right) \mathbf{d}_k

En dos dimensiones es un sistema $2\times2$ . Desarrollado queda así.

\begin{pmatrix} \sum w\, d_x^2 & \sum w\, d_x d_y \\[2pt] \sum w\, d_x d_y & \sum w\, d_y^2 \end{pmatrix} \begin{pmatrix} \partial_x\phi \\[2pt] \partial_y\phi \end{pmatrix} = \begin{pmatrix} \sum w\, \Delta\phi\, d_x \\[2pt] \sum w\, \Delta\phi\, d_y \end{pmatrix}

$\Delta\phi = \phi_k - \phi_C$ . La matriz de la izquierda queda fijada solo por la disposición de los vecinos. Se puede calcular y guardar una sola vez por celda. Las ventajas son grandes. No hay suposición de interpolación en la cara. Es insensible a la distorsión. Con tres o más vecinos siempre tiene solución.

Ponderar por distancia#

El peso $w_k$ determina la exactitud. Una elección común es el inverso de la distancia.

w_k = \frac{1}{\lvert \mathbf{d}_k \rvert^{\,n}}

Con $n=0$ todos los vecinos pesan igual. Un vecino lejano tiene una diferencia de valor grande, y con ella también un error de curvatura grande. Ese vecino tira del ajuste del plano. Al subir $n$ a 1, el vecino lejano pierde voz. El peso recae en los vecinos cercanos, donde la suposición lineal local se cumple mejor. En problemas no viscosos suele usarse $n=1$ .

En la siguiente simulación se pueden manipular los parámetros de forma directa.

angle error 16.8° · |∇φ| 0.64 vs 1.08 true — the reconstructed gradient is skewed

Stencil anisotropy: 0.35Field curvature κ: 1.40Weight exponent n: 0.0Neighbor clustering: 0.60

Al bajar anisotropy cerca de 0.2 para achatar el stencil y subir la curvatura κ, el gradiente reconstruido en naranja se separa del valor real en verde. Aquí conviene subir el weight exponent n de 0 a 1. La flecha naranja vuelve hacia el valor real. Es el resultado de amortiguar el vecino lejano y reducir el error de curvatura.

Python — campo lineal exacto, superficie curva con ponderación#

Se verifica la potencia de los mínimos cuadrados en dos partes. Primero se reconstruye el gradiente del campo lineal $\phi = 1.3x - 0.7y + 2$ sobre un stencil disperso al azar. En teoría el error debe ser cero. Luego, sobre la superficie curva $\phi = e^{0.7x}\cos(0.9y)$ , se compara el caso sin ponderación contra la ponderación $n=1$ .

import numpy as np
 
def weighted_lsq_gradient(offsets, dphi, n):
    # offsets: (K,2) desplazamientos de vecinos d_k; dphi: (K,) diferencias phi_k - phi_C
    w = 1.0 / np.maximum(np.linalg.norm(offsets, axis=1), 1e-12) ** n
    M = (offsets[:, :, None] * offsets[:, None, :] * w[:, None, None]).sum(axis=0)
    b = (offsets * (dphi * w)[:, None]).sum(axis=0)
    return np.linalg.solve(M, b)
 
xc = np.array([0.4, -0.2])
 
# (A) campo lineal — con cualquier stencil, los minimos cuadrados son exactos
lin = lambda x, y: 1.30 * x - 0.70 * y + 2.0
scatter = np.array([[0.45, 0.30], [-0.40, 0.12], [0.05, -0.48], [-0.33, -0.29],
                    [0.90, 0.15], [-0.11, 0.44], [0.28, -0.09]])   # stencil irregular
dlin = np.array([lin(*(xc + d)) - lin(*xc) for d in scatter])
gA = weighted_lsq_gradient(scatter, dlin, 1.0)
print(f"linear field   LSQ = ({gA[0]:+.6f}, {gA[1]:+.6f})   exact = (+1.300000, -0.700000)")
 
# (B) superficie curva — la ponderacion por distancia recorta el error de reconstruccion
phi  = lambda x, y: np.exp(0.7*x) * np.cos(0.9*y)
grad = lambda x, y: np.array([0.7*np.exp(0.7*x)*np.cos(0.9*y),
                             -0.9*np.exp(0.7*x)*np.sin(0.9*y)])
near = np.array([[0.10, 0.05], [-0.08, 0.09], [0.06, -0.10], [-0.09, -0.06], [0.11, 0.10]])
offs = np.vstack([near, [0.90, 0.70]])          # un vecino alejado
dphi = np.array([phi(*(xc + d)) - phi(*xc) for d in offs])
gt   = grad(*xc)
for n in (0.0, 1.0):
    g = weighted_lsq_gradient(offs, dphi, n)
    err = np.linalg.norm(g - gt)
    print(f"curved field   n={n:.0f}  LSQ = ({g[0]:+.3f}, {g[1]:+.3f})   |error| = {err:.4f}")
print(f"curved field   true = ({gt[0]:+.3f}, {gt[1]:+.3f})")

La salida es esta.

linear field   LSQ = (+1.300000, -0.700000)   exact = (+1.300000, -0.700000)
curved field   n=0  LSQ = (+0.885, +0.200)   |error| = 0.0293
curved field   n=1  LSQ = (+0.891, +0.205)   |error| = 0.0222
curved field   true = (+0.911, +0.213)

En el campo lineal, por más torcido que esté el stencil, la exactitud llega hasta el sexto decimal. Esa es la propiedad clave de los mínimos cuadrados. En la superficie curva, un solo vecino alejado agranda el error. La ponderación por distancia con $n=1$ amortigua ese vecino y reduce el error en cerca de una cuarta parte.

Trampas que no hay que pasar por alto al programar#

Que los mínimos cuadrados resistan la distorsión no los hace infalibles. La matriz de la izquierda $M = \sum w_k \mathbf{d}_k \mathbf{d}_k^{\mathsf T}$ puede enfermarse. Si los vecinos se alinean sobre una misma recta, $M$ se acerca a ser singular. El gradiente queda bien definido en la dirección del stencil, pero en la dirección transversal queda casi libre. El número de condición explota.

En la siguiente simulación se puede manipular de forma directa.

cond(M) = 1.4 (λ_max 3.14 / λ_min 2.271) — well-conditioned

Collinearity: 0.20Neighbor count: 6

Al subir Collinearity hasta cerca de 0.9, los vecinos se agrupan en una sola fila y cond(M) se dispara de decenas hasta el infinito. La dirección en que la elipse se estira larga y delgada es justo la dirección en que el gradiente se vuelve inexacto. Es una situación que ocurre a menudo en las celdas de capa límite pegadas a la pared.

Tres cosas a cuidar en la práctica. Primero, en las celdas de frontera los vecinos se cargan hacia un lado. Conviene agregar al stencil un vecino del lado opuesto a la pared o el centro de una cara para bajar el número de condición. Segundo, no reconstruir $M$ en cada paso: si la malla es fija, se guarda de antemano la inversa (o la descomposición QR). Tercero, ajustar la ponderación por distancia $n$ según el problema. Si el término viscoso domina, a veces $n=2$ resulta mejor.

Queda en una línea. Green–Gauss es barato pero débil ante la distorsión, y los mínimos cuadrados ponderados resisten la distorsión pero se enferman si el stencil se alinea en una sola fila — hay que mirar primero la malla y luego elegir el método.