Los extremos falsos que una reconstrucción de segundo orden inventa en mallas no estructuradas — Limitadores Barth–Jespersen y Venkatakrishnan

El solver venía bien. Al activar el segundo orden, la densidad cayó a negativo. Delante del choque, el valor de cara reconstruido quedaba por encima del promedio de celda. El upwind de primer orden, en la misma malla, seguía sano. Subir el orden fue lo que mató la corrida. Una recta dentro de una celda no estructurada inventa extremos que ningún vecino tiene. Barth–Jespersen corta esos extremos. Venkatakrishnan redondea el corte.

La reconstrucción inventa extremos que nunca estuvieron ahí#

En el método de volúmenes finitos (FVM, una discretización que lleva promedios de celda como incógnitas), el segundo orden es una recta dentro de la celda.

$u_i$ es el promedio en la celda $i$. $\nabla u_i$ es el gradiente reconstruido, $\mathbf{x}_f$ es el centro de la cara y $\Phi_i\in[0,1]$ es el limitador.

Con $\Phi_i=1$ el limitador está apagado. Mínimos cuadrados o Green–Gauss arman el gradiente promediando vecinos. Junto a un salto, ese promedio queda demasiado empinado. Al extender la recta hasta la cara aparece un máximo o un mínimo nuevo, uno que ningún vecino guardaba.

Al subir la ganancia del gradiente, la reconstrucción sin limitar (Φ=1) atraviesa la banda gris — el máximo y el mínimo de los vecinos — y pasa a rojo.

Mientras la línea permanece dentro de la banda, la Φ de esa celda vale 1. Solo las celdas junto al salto se acortan primero. Con la misma ganancia, al pasar el limitador a Barth–Jespersen la línea deja de salir de la banda.

El extremo falso contamina el flujo del paso siguiente. La oscilación termina en densidad negativa y energía negativa. Un valor de cara reconstruido no debe salir del rango que sus vecinos ya tienen. Ese es el principio de máximo local (LMP). $\Phi_i$ es el escalar que recorta el gradiente para obedecerlo.

Barth–Jespersen: enjaulado por el máximo y el mínimo de los vecinos#

Barth y Jespersen (1989) fijan primero el rango permitido. Es el máximo y el mínimo de la celda $i$ y de sus vecinos.

$N(i)$ es el conjunto de vecinos por cara de la celda $i$. En cada cara se lee el incremento sin limitar $u_f-u_i$ y se forma un factor por cara $\phi_f$.

El limitador de la celda es el valor más conservador entre todas las caras.

Dentro del rango, $\phi_f=1$ y el segundo orden se queda. Una cara que se pasaría conserva solo la pendiente justa para tocar el borde. Cerca de un extremo $\Phi_i\to 0$, y el esquema retrocede al upwind de primer orden.

Para un incremento positivo el factor es $\phi_f=\min(1,(u_i^{\max}-u_i)/(u_f-u_i))$. En el instante en que el numerador es el menor, el valor de cara se pega a $u_i^{\max}$ y se detiene. Esa razón explica por qué la línea roja de la figura de arriba se quiebra en el borde gris.

La segunda trampa: la no diferenciabilidad#

Barth–Jespersen sostiene la monotonía. Al meterlo en un solver estacionario, el residuo se estanca y vibra en cierto nivel. La causa es el $\min$ y la división. $\phi_f$ no es diferenciable respecto a la solución. Qué cara aporta el mínimo cambia de una iteración a otra. El limitador se enciende y se apaga, y el residuo rebota. El jacobiano de un solver implícito detesta ese quiebre.

En regiones suaves, Φ basta con valer 1. Solo debe activarse en una discontinuidad real, y el paso entre encendido y apagado tiene que ser suave. La expresión siguiente redondea el quiebre agudo.

Venkatakrishnan: recortando suavemente#

Venkatakrishnan (1993) reemplazó el quiebre de $\min(1,y)$ por una función racional. El limitador Michalak–Ollivier-Gooch que menciona la fuente pertenece a la misma familia. El factor por cara se escribe

$\Delta_-=u_f-u_i$ es el incremento sin limitar. $\Delta_+$ es $u_i^{\max}-u_i$ o $u_i^{\min}-u_i$, según el signo. La cantidad que se elige a mano es $\epsilon^2$.

$K$ es el parámetro de ajuste y $\Delta x$ es el tamaño de malla. Cuando la variación es menor que $\epsilon$, $\phi_f\to 1$ y el limitador se apaga. Al subir $K$ el umbral sube, así que el limitador se suelta sobre una región más amplia. Demasiado alto, y las oscilaciones junto al choque se escapan. Con $K$ en 0, la función vuelve a Barth–Jespersen.

El código toma el tamaño de celda como 1 y calcula $\epsilon^2=K^3$. Con el valor base $K=0.3$, $\epsilon^2=0.027$. Si $\epsilon=0$ y $\Delta_+=\Delta_-=\Delta$, justo donde la reconstrucción cae exactamente sobre el límite permitido, la función racional vale $(3\Delta^2)/(4\Delta^2)=3/4$. Barth–Jespersen devolvería $\phi_f=1$ en ese mismo punto. Venkatakrishnan tira un cuarto de la pendiente antes de que la cara toque el límite. Por eso la esquina se embota un paso antes.

El intervalo que se prueba por costumbre va de $K\approx 0.1$ a $K\approx 5$. Al meter esos extremos en el $\epsilon^2=K^3$ del código salen 0.001 y 125. La razón entre umbrales es 125000. Un movimiento chico de $K$ va al cubo, así que el conjunto de celdas excusadas por suaves crece rápido. Bajar $K$ de 0.1 deja $\epsilon^2$ prácticamente en 0, y regresa el quiebre de Barth–Jespersen.

En la figura de arriba, al activar Venkatakrishnan y mover $K$, Φ se desliza sin escalones.

Python — carrera de tres limitadores en advección escalar#

Dominio periódico. Una onda cuadrada y una campana gaussiana se transportan juntas (advección). Sin limitar (Fromm), Barth–Jespersen y Venkatakrishnan corren en paralelo. Importa el mínimo final. Por debajo de cero, se cavó un valle falso de esa profundidad.

La pendiente de celda es $s=(u_{i+1}-u_{i-1})/2$. La cara recibe el incremento $0.5s$, de modo que una diferencia de 1 entre vecinos da $s=1/2$ y un incremento de cara de 0.25. Si ese 0.25 supera el margen permitido, Φ baja de 1. Una diferencia central unidimensional solo ve a los vecinos izquierdo y derecho. En una celda no estructurada, los pesos de ese mismo gradiente son otra decisión: la que los mínimos cuadrados móviles reajustan celda por celda. El limitador se engancha después de que el gradiente ya existe.

import numpy as np
 
NX, A, CFL = 200, 1.0, 0.4
dx = 1.0 / NX
dt = CFL * dx / A
 
def init_profile():
    x = (np.arange(NX) + 0.5) / NX
    u = np.where((x > 0.1) & (x < 0.3), 1.0, 0.0)   # onda cuadrada
    u += np.exp(-((x - 0.65) / 0.06) ** 2)          # campana gaussiana
    return u
 
def cell_slope(u):
    return (np.roll(u, -1) - np.roll(u, 1)) / 2.0    # pendiente central (Fromm)
 
def barth_jespersen_phi(u, s):
    up, um = np.roll(u, -1), np.roll(u, 1)
    umax = np.maximum(u, np.maximum(up, um))
    umin = np.minimum(u, np.minimum(up, um))
    phi = np.ones_like(u)
    for du in (0.5 * s, -0.5 * s):                   # dos caras
        f = np.ones_like(u)
        pos, neg = du > 1e-12, du < -1e-12
        f[pos] = np.minimum(1.0, (umax[pos] - u[pos]) / du[pos])
        f[neg] = np.minimum(1.0, (umin[neg] - u[neg]) / du[neg])
        phi = np.minimum(phi, f)
    return np.clip(phi, 0.0, 1.0)
 
def venkatakrishnan_phi(u, s, K=0.3):
    up, um = np.roll(u, -1), np.roll(u, 1)
    umax = np.maximum(u, np.maximum(up, um))
    umin = np.minimum(u, np.minimum(up, um))
    eps2 = K ** 3                                    # (K*h)^3, tamaño de celda h=1
    phi = np.ones_like(u)
    for du in (0.5 * s, -0.5 * s):
        d = np.where(du > 0, umax - u, umin - u)
        num = (d * d + eps2) * du + 2 * du * du * d
        den = d * d + 2 * du * du + d * du + eps2
        f = np.where(np.abs(du) < 1e-12, 1.0, num / (du * den))
        phi = np.minimum(phi, f)
    return np.clip(phi, 0.0, 1.0)
 
def muscl_rhs(u, limiter):
    s = cell_slope(u)
    phi = limiter(u, s) if limiter else np.ones_like(u)
    uL = u + 0.5 * phi * s              # estado izquierdo en la cara i+1/2
    flux = A * uL                       # flujo upwind (a > 0)
    return -(flux - np.roll(flux, 1)) / dx
 
def advance_muscl(u, limiter):         # SSP-RK2
    k1 = muscl_rhs(u, limiter)
    u1 = u + dt * k1
    k2 = muscl_rhs(u1, limiter)
    return 0.5 * (u + u1 + dt * k2)
 
def run_limiter_race(steps=160):
    fields = {"none": init_profile(), "bj": init_profile(), "venk": init_profile()}
    lims = {"none": None, "bj": barth_jespersen_phi, "venk": venkatakrishnan_phi}
    for _ in range(steps):
        for k in fields:
            fields[k] = advance_muscl(fields[k], lims[k])
    for k, u in fields.items():
        print(f"{k:5s}  min={u.min():+.4f}  max={u.max():.4f}")
 
run_limiter_race()

La corrida imprime esto.

none   min=-0.0417  max=1.0231
bj     min=+0.0000  max=1.0000
venk   min=-0.0004  max=1.0009

El mínimo sin limitar, -0.0417, es lo que una onda cuadrada de altura 1 se hundió bajo cero. El máximo, 1.0231, es lo que se pasó de esa misma altura. Barth–Jespersen clava mínimo y máximo en $[0,1]$. El mínimo de Venkatakrishnan, -0.0004, es cerca de 1/104 de ese valle. Lo que supera a 1 baja de 0.0231 a 0.0009.

La distancia en 160 pasos es $A \cdot dt \cdot 160$. Con $dx=1/200$ y CFL 0.4, $dt=0.002$, así que la distancia es 0.32. Son 64 celdas de 200. El intervalo de la onda cuadrada $(0.1, 0.3)$ pasa a $(0.42, 0.62)$, y el centro gaussiano 0.65 llega a 0.97. Las dos estructuras se mueven a la misma velocidad, de modo que la posición relativa se mantiene. El mínimo y el máximo impresos no son un choque entre ondas. Cada esquina produjo su propio número. La longitud 0.2 de la onda cuadrada son 40 celdas. La escala 0.06 de la gaussiana son 12 celdas. Unas pocas celdas junto a esas esquinas fijan el mínimo y el máximo.

Basta el borde de salida de la onda cuadrada. El esquema sin limitar (rojo) hace crecer rizos por debajo de cero. Los dos limitadores se mantienen monótonos. Venkatakrishnan (ámbar) queda un poco más romo en la esquina que Barth–Jespersen (cian). Apunta en la misma dirección que el álgebra, que con $\epsilon=0$ conservaba solo 3/4 de la pendiente. El $K=0.3$ del código no borra el quiebre. Solo reduce el exceso a un nivel de 0.0009.

Al encender el limitador en la práctica#

Cómo se agrupan los vecinos fija el rango permitido. Tomar el $\max/\min$ de Barth–Jespersen sobre vecinos por cara o sobre vecinos por vértice cambia el resultado. En mallas no estructuradas, los vecinos por vértice — toda celda que comparte un vértice — reducen el sesgo direccional y equilibran el limitador.

Después viene $K$. Demasiado pequeño, y la convergencia se estanca como con Barth–Jespersen. Demasiado grande, y las oscilaciones del choque se cuelan. Como $\epsilon^2$ lleva $\Delta x^3$, refinar la malla hace que el limitador se dispare más seguido con el mismo $K$.

En un sistema vectorial (las ecuaciones de Euler), limitar cada variable conservada por separado puede crear oscilaciones nuevas por desajuste entre componentes. Proyectar primero a variables características (characteristic variable) y limitar después es más seguro. Cuesta más.

Si el residuo se detiene, el gradiente no es el primer sospechoso#

Cuando un solver de segundo orden muere en densidad negativa, lo primero es ver si el valor de cara salió del rango de los vecinos. Si salió, el problema está en el conjunto de vecinos que armó Φ, no en la rutina del gradiente. Los vecinos por cara solos pierden extremos en la diagonal. La reconstrucción vuelve a filtrarse por esa rendija.

Si en un cálculo estacionario el residuo solo da pasos en el sitio, se abre $K$. $K=0$ es el quiebre de Barth–Jespersen, así que el jacobiano salta en cada iteración. Donde un salto lleva Φ a 0, el orden de esa celda se queda en 1. ¿La región suave sigue cerca del segundo orden? Si sigue, $\epsilon^2$ sí dejó libre esa región. Qué números aparecen cuando una sola corrida mide las dos cosas a la vez sigue abierto en el registro que pone un choque y un vórtice suave en la misma prueba de orden.