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

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$. Las dos estimaciones son Green–Gauss y mínimos cuadrados. La figura de abajo muestra dónde la distorsión tumba a una, y cuánto devuelve el peso por distancia.

Al bajar anisotropy cerca de 0.2 el stencil se aplana. Al subir el weight exponent de 0 a 1, la flecha naranja gira de nuevo hacia el valor real verde.

Si curvature κ sube a la vez, la separación entre las flechas se abre. En un stencil plano, lo que importa es si el peso recupera la dirección.

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.

$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.

$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.

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

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

$\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.

La interpolación en la cara se rompe en el mismo sitio que Green–Gauss frente a mínimos cuadrados en una celda sesgada. Con tres vecinos las ecuaciones normales tienen solución, pero si esos tres quedan en una recta la respuesta transversal sale vacía.

Ponderar por distancia#

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

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$. Ese valor recorta el voto del vecino lejano y baja 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 un campo lineal el stencil puede estar todo lo torcido que se quiera y el gradiente sigue bien hasta el sexto decimal. La salida (+1.300000, -0.700000) es ese valor. Una función lineal hace coincidir el producto punto del desplazamiento y el gradiente en cada vecino, así que el punto irregular (0.90, 0.15) no aporta error. La distorsión no empuja la solución en este campo.

En una superficie curva sí. El gradiente real es (+0.911, +0.213). Con n=0, la reconstrucción (+0.885, +0.200) se queda corta en 0.026 en x y 0.013 en y. Con n=1, (+0.891, +0.205) se queda corta en 0.020 y 0.008. La norma baja de 0.0293 a 0.0222. La diferencia, 0.0071, es cerca del 24 por ciento, el recorte de un cuarto. El vecino lejano (0.90, 0.70) tiene componentes mucho mayores que la escala 0.1 de los cinco puntos cercanos. Con n=0 ese punto vota igual que los demás y arrastra el plano. Con n=1 la distancia le recorta el voto.

Si los vecinos se alinean, la inversa no sirve#

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.

Al subir Collinearity hasta cerca de 0.9, los vecinos se juntan en una sola fila y cond(M) salta de decenas al infinito.

El eje largo y delgado de la elipse es el eje que el gradiente no logra fijar. Las celdas de capa límite pegadas a la pared tienen esa forma. La diagonal es una suma de cuadrados de los desplazamientos, así que un stencil estirado en una sola dirección deja la otra diagonal vacía. En una celda de frontera los vecinos se cargan del lado del fluido. Agregar un vecino del lado opuesto a la pared, o el centro de una cara, rellena un poco ese eje vacío y el número de condición baja. Si la malla está fija, $M$ no cambia aunque cambie el campo de velocidad. Conviene guardar la inversa, o la descomposición QR, en lugar de reconstruirla en cada paso. El exponente $n$ depende del problema. Cuando el término viscoso es grande, a veces $n=2$ resulta mejor.

Los mínimos cuadrados móviles rehacen ese peso vecino por vecino. Green–Gauss es barato y débil ante el sesgo. Los mínimos cuadrados ponderados aguantan el sesgo, hasta que el stencil se alinea y la matriz enferma primero. Antes de elegir el método, hay que ver si los vecinos de esa celda cubren el plano.