[Reseña] Cuando WENO difumina los pequeños vórtices — Fu (2019) esquema TENO de baja disipación en volúmenes finitos

Con ese mismo 10 por ciento, los pesos normalizados de TENO se quedan cerca de un cuarto y el estencil ancho sigue activo.

En Stanford, Lin Fu corría una DNS de turbulencia y vio algo raro. WENO-JS de quinto orden atrapaba los choques. También emborronaba vórtices pequeños lejos de esos choques. El esquema decía quinto orden. La resolución en la mano se parecía más a una diferencia central de tercer orden.

Al mover el deslizador de $C_T$ de $10^{-5}$ a $10^{-3}$ o a $10^{-7}$ , el corte se pone exigente o se afloja, y el mismo bulto se recorta distinto en naranja y en cian.

N cells: 120

CFL: 0.40

TENO cutoff C_T: 1e-5

WENO-JSTENO-5L2 err — WENO: 0.0000 | TENO: 0.0000

Tras un período, conviene mirar si los dos bordes del pulso cuadrado se emborronan solo una celda o dos. También si a $10^{-7}$ aparece un anillo pequeño junto al pulso cian.

Información del artículo y contexto#

Autor: Lin Fu (Stanford Turbulence Research Center)
Revista: Computer Physics Communications 235 (2019) 25–39
DOI: 10.1016/j.cpc.2018.10.009
Palabras: TENO, High-order schemes, Shock-capturing, Low-dissipation, Finite-volume method

El costo oculto de WENO-JS#

El WENO-JS clásico (Jiang–Shu, 1996) mide un indicador de suavidad $\beta_k$ en cada estencil pequeño $S_0, S_1, S_2$ y los mezcla con pesos no lineales.

w_k = \frac{\alpha_k}{\sum_j \alpha_j}, \qquad \alpha_k = \frac{d_k}{(\beta_k + \varepsilon)^2}

$d_k$ es el peso lineal óptimo que devuelve el quinto orden. $\varepsilon$ , cerca de $10^{-6}$ , mantiene el denominador lejos de cero. Esos pesos son continuos. Un temblor mínimo del estencil mueve $w_k$ en una zona que debía quedarse en $w_k=d_k$ , y el peso se fuga hacia un estencil de orden menor.

Una onda suave pierde amplitud y acumula error de fase en decenas de pasos. En LES, donde choques y vórtices turbulentos comparten el dominio, esa disipación es fatal.

Dos ideas detrás de TENO: separación de escalas fuerte y corte nítido#

Fu cambió dos cosas.

Cambió el conjunto de estenciles. Junto a los tres estenciles de 3 puntos $S_0, S_1, S_2$ hay un estencil grande de 5 puntos $S_3$ . Si $S_3$ está limpio, se detiene. Los otros no se usan.

Los pesos se volvieron binarios. La prueba de suavidad es el corte $C_T$ .

\gamma_k = \left(1 + \frac{\tau_K}{\beta_k + \varepsilon}\right)^q, \qquad \chi_k = \frac{\gamma_k}{\sum_j \gamma_j}

\delta_k = \begin{cases} 0, & \chi_k < C_T \\ 1, & \text{otherwise} \end{cases}

$\tau_K=|\beta_0-\beta_2|$ es la referencia global. $q=6$ estira el cociente. $\delta_k$ es 0 o 1. Un estencil se descarta del todo o se queda en su peso óptimo. No hay un fundido continuo. Las zonas suaves conservan el esquema lineal de quinto orden sobre $S_3$ . Solo cerca de una discontinuidad el esquema elige estenciles al modo ENO. Si un solo corte puede cargar con toda onda es otra pregunta. El argumento de que cada familia de ondas necesita su propia reconstrucción duda del intento.

De pesos continuos a un veredicto 0/1#

Si $S_3$ es suave, $\delta_3=1$ , la reconstrucción es la fórmula lineal de quinto orden.

u_{i+1/2}^L = \frac{1}{60}\bigl(2\bar u_{i-2} - 13\bar u_{i-1} + 47\bar u_i + 27\bar u_{i+1} - 3\bar u_{i+2}\bigr)

$\bar u_j$ es el promedio de celda. La fórmula es del todo lineal. Cuando $C_T$ está cerca de cero, disipa casi como una diferencia central. Si $\delta_3=0$ , solo se mezclan los estenciles pequeños que sobrevivieron, con $d_k$ .

u_{i+1/2}^L = \sum_{k=0}^{2} w_k\, u_{k,i+1/2}^L, \qquad w_k = \frac{d_k\, \delta_k}{\sum_j d_j\, \delta_j}

Esa es la propiedad ENO. Un estencil que cruza una discontinuidad tiene $\delta_k$ apagado, y el esquema pasa a otro estencil.

Reconstrucción de estencil de 5° orden en NumPy#

El código es corto. Con cinco puntos arma tres indicadores de suavidad y tres polinomios candidatos. TENO corta con $C_T$ . WENO-JS, la comparación, mezcla por inversos de cuadrados. Esta función es la rama TENO.

import numpy as np
 
def teno_beta_indicators(um2, um1, u0, up1, up2):
    # Indicadores de suavidad de Jiang-Shu para tres estenciles de 3 puntos
    b0 = (13/12)*(um2 - 2*um1 + u0)**2 + 0.25*(um2 - 4*um1 + 3*u0)**2
    b1 = (13/12)*(um1 - 2*u0 + up1)**2 + 0.25*(um1 - up1)**2
    b2 = (13/12)*(u0 - 2*up1 + up2)**2 + 0.25*(3*u0 - 4*up1 + up2)**2
    return b0, b1, b2
 
def teno_weighted_reconstruction(u, i, CT=1e-5, q=6):
    N = len(u)
    um2, um1, u0 = u[(i-2) % N], u[(i-1) % N], u[i]
    up1, up2 = u[(i+1) % N], u[(i+2) % N]
    b0, b1, b2 = teno_beta_indicators(um2, um1, u0, up1, up2)
    b3 = max(b0, b1, b2)
    tau = abs(b0 - b2)
    eps = 1e-40
    g = np.array([(1 + tau/(b + eps))**q for b in (b0, b1, b2, b3)])
    chi = g / g.sum()
    delta = (chi >= CT).astype(float)
    if delta[3] == 1:  # S3 suave -> usar 5° orden lineal tal cual
        return (2*um2 - 13*um1 + 47*u0 + 27*up1 - 3*up2) / 60
    d = np.array([1/10, 6/10, 3/10])
    alpha = d * delta[:3]
    if alpha.sum() < 1e-14:
        return u0  # celda donadora como último recurso
    w = alpha / alpha.sum()
    p0 = (2*um2 - 7*um1 + 11*u0) / 6
    p1 = (-um1 + 5*u0 + 2*up1) / 6
    p2 = (2*u0 + 5*up1 - up2) / 6
    return w[0]*p0 + w[1]*p1 + w[2]*p2

Con un paso de Euler y fronteras periódicas, un bulto gaussiano suave vuelve casi intacto con TENO y cerca de un 30 por ciento más bajo con WENO-JS. Mismo esqueleto. Solo cambian los pesos.

Se sustituye un temblor pequeño en los pesos. Los ideales son un décimo, seis décimos y tres décimos. Si los tres β valen 1, WENO vuelve ahí. Subir solo el tercero un 10 por ciento pone β en 1, 1, 1.1. El tercer peso sin normalizar es 0.3 sobre 1.21, cerca de 0.248. Normalizados: cerca de 0.105, 0.633 y 0.262. La brecha L1 es cerca de 0.077.

Los mismos números entran en TENO. Tau es 0.1 y q es 6. El indicador del estencil grande es el máximo de los tres β, y se normalizan cuatro juntos. Chi queda cerca de 0.256, 0.256, 0.244 y 0.244. Frente a $10^{-5}$ eso está cuatro órdenes arriba. Todo veredicto se queda en 1 y se usa la fórmula lineal de quinto orden.

Épsilon en la fórmula anda cerca de $10^{-6}$ . Dentro de la función es $10^{-40}$ . Cerca de β igual a 1 ninguno aparece. Solo un estencil con β en 0 parte el denominador. Los coeficientes 2, −13, 47, 27 y −3 suman 60, así que un campo constante no cambia.

Qué queda cuando C_T se mueve#

La advección de arriba empuja una sola condición inicial, un bulto de seno al cuadrado más un pulso cuadrado. El naranja es WENO-JS. El cian es TENO-5. Con $C_T=10^{-5}$ , tras un período, los dos bordes del pulso se emborronan una celda o dos. Solo el bulto seno naranja baja de un modo visible. Subir $C_T$ a $10^{-3}$ vuelve exigente el corte, los estenciles pequeños ganan a menudo junto al pulso, y TENO disipa más. Bajarlo a $10^{-7}$ hace que TENO insista en el esquema lineal de quinto orden, y aparece un anillo pequeño junto al pulso. $C_T$ es la balanza entre estabilidad y baja disipación.

Notas críticas — La trampa del ajuste#

Unos sitios que el artículo no afirma en voz alta.

$C_T$ se pega al problema. Fu recomienda $10^{-5}$ . Con un choque fuerte y una onda acústica muy débil, ese valor levanta un poco la onda acústica, y $10^{-6}$ hace sonar el entorno del choque. Los grupos que lo portaron a OpenFOAM o SU2 escriben que reajustan el valor en cada caso.

En volumen finito, la suavidad $\beta_3$ de $S_3$ sale cara. Las derivadas del polinomio de 5 puntos llenan unas seis páginas. La forma en diferencias finitas es mucho más corta. Conservar el cambio de flujo de Riemann de baja disipación de este artículo pide el arreglo de volumen finito.

El encaje con la integración en el tiempo está vacío. SSP-RK3 va bien, tapado por la exactitud espacial. Mezclado con SSP-RK2 o un método implícito, se reporta que el salto de $\delta_k$ molesta la convergencia de Newton.

En la práctica, el estudio de malla se corre con un $C_T$ fijo. Conviene revisar si el orden espacial vuelve de verdad a 5. Si no vuelve, $C_T$ probablemente es tan grande que TENO se queda todo el tiempo en el modo no lineal.

La verificación que se rompe si C_T queda fijo#

Al refinar la malla con $C_T$ congelado en $10^{-5}$ , el quinto orden desaparece en el momento en que chi de una celda suave cae bajo el corte. El temblor del 10 por ciento de arriba no cruzó el corte. Un choque junto a una onda muy débil es otra historia.

Si entonces una prueba de orden vuelve por debajo de 1 en lugar de 5, el corte se comió el orden, no el operador espacial. La verificación en la que un mismo código reportó orden 0.79 y orden 2.00 a la vez pregunta qué está midiendo de verdad el caso. ¿En cuál de esos dos números cae TENO con un $C_T$ congelado?