[Reseña] Dibujar choques con tanh — Implementando la reconstrucción TENO-THINC

Subir el orden formal a quinto, sexto u octavo no impide que el contacto se difumine con el tiempo. WENO o TENO, igual. Dibujar un salto con un polinomio hace eso. El artículo de 2023 de Takagi y colegas se aparta de ese camino. El TENO de alto orden se queda con las regiones suaves. Solo las celdas discontinuas cambian a una tanh (THINC).

Al cambiar la reconstrucción entre Polynomial, THINC e Hybrid, los mismos promedios de celda se parten en una curva S, un salto vertical o una mezcla.

Polynomial se pasa un poco por arriba y por abajo del salto. THINC queda casi vertical dentro de una celda. Hybrid reserva la tanh para la celda del salto.

El artículo en una página#

• Título: High-Order Low-Dissipation Shock-Resolving TENO-THINC Schemes for Hyperbolic Conservation Laws
• Autores: Shinichi Takagi, Hiro Wakimura, Lin Fu, Feng Xiao
• Publicación: Communications in Computational Physics, 2023
• En una línea: Acoplar el TENO de puntos pares (seis y ocho puntos) con la reconstrucción THINC para lograr baja disipación en las regiones suaves y resolución sub-celda en las discontinuidades al mismo tiempo.

Lo que sale gratis es decidir cuándo cambiar. TENO ya marca cada plantilla candidata como suave o no mediante $\delta_k \in \{0,1\}$. Reutilizar esas banderas indica si una celda esconde una discontinuidad, y THINC se activa solo ahí. El único parámetro nuevo es la pendiente $\beta$ de THINC.

Por qué incluso los polinomios de alto orden difuminan el contacto#

La reconstrucción polinómica dibuja una curva suave por naturaleza. Cuando una curva debe pasar por los puntos de promedio de celda a ambos lados de un salto casi vertical, enfrenta una disyuntiva: permitir oscilación (sobreimpulso), o matar la pendiente y difuminar.

TENO y WENO eligen suprimir la oscilación. Descartan las plantillas que cruzan la discontinuidad, imponiendo la propiedad ENO. Eso compra estabilidad, pero un poco de disipación numérica se acumula en cada paso. Una discontinuidad de contacto (una superficie donde la velocidad y la presión coinciden pero la densidad salta) no se auto-afila por un mecanismo de presión. Así que tras una advección prolongada su ancho se desangra a lo largo de decenas de celdas. Subir el orden apenas reduce este difuminado.

Las tres reconstrucciones de arriba son esa elección. Mismos promedios, curvas distintas. Donde TENO mata menos los vórtices pequeños está en TENO de baja disipación. Un contacto que se derrama por decenas de celdas no se detiene solo porque baje la disipación.

THINC — la tanh dibuja un salto sub-celda#

THINC (Tangent of Hyperbola for INterface Capturing) reemplaza el polinomio por una tanh monótona. La reconstrucción sobre la coordenada local $X \in [-0.5, 0.5]$ en la celda $i$ es

donde $f_a = (f_{i+1}+f_{i-1})/2$ es la media de vecinos, $f_d = (f_{i+1}-f_{i-1})/2$ es la mitad del salto, $\beta$ es la pendiente de la tanh y $d_i$ es la posición del centro del salto.

$d_i$ no es una variable libre. La restricción de conservación del promedio de celda (la reconstrucción debe integrar de vuelta a $f_i$) lo fija. La solución es analítica:

con $T_1 = \tanh(\beta/2)$, $T_2 = \tanh(\alpha\beta/2)$, $\alpha = (f_i - f_a)/f_d$. Aquí $\alpha$ es la posición normalizada de la media de celda entre sus vecinos. Si la celda es no monótona ($(f_{i+1}-f_i)(f_i-f_{i-1}) \le 0$), la tanh carece de sentido, así que retrocede a primer orden ($\hat f_{i+1/2} = f_i$).

β: la perilla que sostiene el filo#

Un solo $\beta$ fija el carácter de THINC: pequeño es suave, grande es empinado. El artículo registra una correspondencia útil.

• $\beta \approx 1.1$ : espectro TVD al nivel del limitador van Leer
• $\beta \approx 1.3$ : nivel Superbee (tendencia sobre-compresiva)
• $\beta = 1.6 \sim 2.0$ : régimen de captura de choques afilada

El artículo adopta $\beta = 1.8$. Subirlo demasiado escalona hasta las curvas suaves (squaring) y aplana de forma artificial. En la figura de arriba, al mover $\beta$ de 1.0 a 2.4 el salto se yergue y la lectura de overshoot se mueve con él.

Con $\beta = 1.8$, $\beta/2 = 0.9$. Un salto centrado en la celda tiene $\alpha = 0$, luego $d_i = 0$. En el borde derecho $X = 0.5$, $\tanh(\beta(X - d_i)) = \tanh(0.9) \approx 0.716$. Solo vuelve cerca del 72% del salto vecino dentro de una celda. Con $\beta = 1.1$, $\tanh(0.55) \approx 0.50$, tan suave como van Leer. Con $\beta = 2.4$, $\tanh(1.2) \approx 0.83$. El valor 1.8 queda entre una pared vertical y el squaring.

$\alpha = (f_i - f_a)/f_d$. Si la media de celda cae a mitad de camino entre los vecinos, $\alpha = 0$ y el centro del salto es 0. Si la media se pega al vecino derecho, $\alpha$ se acerca a 1, $(1 - T_2/T_1)$ dentro del logaritmo se encoge y $d_i$ se desplaza en positivo. Si $|\alpha| \ge 1$ o $(f_{i+1} - f_i)(f_i - f_{i-1}) \le 0$, el código suelta la tanh y retrocede a $\hat f = f_i$. Las celdas de pico se quedan con el polinomio.

El top-hat inicial está entre $0.35$ y $0.65$. Un ancho de $0.30$ con $N = 100$ son 30 celdas. $dx = 0.01$ y $\Delta t = 0.45 \times dx = 0.0045$. Cinco vueltas son $5 / 0.0045 \approx 1111$ pasos, con CFL 0.45. Después, el ancho de minmod es 23 celdas y el pico 0.84. Siete celdas caen por debajo de 0.5 y la altura pierde 0.16. THINC-BVD mantiene pico 1.00 y ancho 30. Las dos líneas de comentario del código son esa cuenta.

TENO lo encuentra, THINC lo rellena — el híbrido BVD#

Aplica THINC en todas partes y hasta las ondas suaves se rompen. Por eso debe activarse solo en celdas discontinuas. El artículo construye un indicador por celda $\zeta_i$ a partir de los $\delta_k$ que el pesado de TENO ya calculó.

$\zeta_i = 1$ señala una discontinuidad a la derecha, $-1$ a la izquierda. Cuando los signos de $\zeta$ de los vecinos se enfrentan a través de la celda $i$, se juzga que esa celda realmente contiene una discontinuidad y se cambia por THINC.

El BVD (Boundary Variation Diminishing) práctico es aún más simple. Calcula ambas reconstrucciones y conserva la que más reduzca la diferencia entre los valores reconstruidos izquierdo y derecho en la interfaz de celda (la variación de frontera). La implementación de abajo sigue esta prueba BVD.

Python — girar un top-hat una vuelta#

Advecta un top-hat bajo advección lineal $u_t + u_x = 0$ con fronteras periódicas. Compara minmod TVD contra THINC-BVD en la misma malla por el mismo tiempo.

import numpy as np
 
def minmod(p, q):
    return np.where(p * q <= 0, 0.0, np.where(np.abs(p) < np.abs(q), p, q))
 
def thinc_right_edge(u, beta=1.8):
    # Ecs. 2.19-2.20: valor reconstruido izquierdo en la cara derecha (i+1/2) de la celda i
    um, up = np.roll(u, 1), np.roll(u, -1)
    fa, fd = 0.5 * (up + um), 0.5 * (up - um)
    mono = ((up - u) * (u - um) > 0) & (np.abs(fd) > 1e-12)
    alpha = np.where(mono, (u - fa) / np.where(mono, fd, 1.0), 0.0)
    mono = mono & (np.abs(alpha) < 1.0)
    T1, T2 = np.tanh(beta / 2), np.tanh(alpha * beta / 2)
    d = np.where(mono, (1 / (2 * beta)) * np.log((1 - T2 / T1) / (1 + T2 / T1)), 0.0)
    vR = fa + fd * np.tanh(beta * (0.5 - d))
    return np.where(mono, vR, u)  # no monótona -> primer orden
 
def reconstruct(u, beta, use_thinc):
    um, up = np.roll(u, 1), np.roll(u, -1)
    s = minmod(u - um, up - u)
    poly_vR = u + 0.5 * s
    if not use_thinc:
        return poly_vR
    poly_vL = u - 0.5 * s
    t = thinc_right_edge(u, beta)
    # BVD: conservar la reconstruccion que reduce el salto en la cara derecha
    jump_poly = np.abs(poly_vR - np.roll(poly_vL, -1))
    jump_thinc = np.abs(t - np.roll(poly_vL, -1))
    return np.where(jump_thinc < jump_poly, t, poly_vR)
 
def rhs(u, dx, beta, use_thinc):
    vR = reconstruct(u, beta, use_thinc)        # a=1, asi que F_{i+1/2} = vR
    return -(vR - np.roll(vR, 1)) / dx
 
def advect_tophat(N=100, revolutions=5.0, beta=1.8, use_thinc=True):
    dx, x = 1.0 / N, (np.arange(N) + 0.5) / N
    dt = 0.45 * dx
    u = ((x > 0.35) & (x < 0.65)).astype(float)   # top-hat inicial
    for _ in range(int(revolutions / dt)):
        u1 = u + dt * rhs(u, dx, beta, use_thinc)            # SSP-RK2
        u = 0.5 * u + 0.5 * (u1 + dt * rhs(u1, dx, beta, use_thinc))
    return x, u
 
x, u_tvd = advect_tophat(use_thinc=False)
_, u_thinc = advect_tophat(use_thinc=True)
print("TVD   peak=%.3f  width(>0.5)=%d" % (u_tvd.max(), (u_tvd > 0.5).sum()))
print("THINC peak=%.3f  width(>0.5)=%d" % (u_thinc.max(), (u_thinc > 0.5).sum()))
# Salida de ejemplo:
# TVD   peak=0.84  width(>0.5)=23
# THINC peak=1.00  width(>0.5)=30

Tras cinco vueltas, minmod TVD se hunde a un pico de 0.84 y se estrecha. THINC-BVD mantiene el pico 1.00 y conserva casi su ancho inicial (30 celdas). Ni siquiera un TENO de mayor orden puede frenar del todo este difuminado, por la naturaleza polinómica misma — ese es el punto de partida del artículo.

Reconstrucción y advección de largo plazo, en vivo#

Al avanzar revolutions, la curva roja de minmod y la cian de THINC recorren el mismo tiempo por separado.

Conforme crece revolutions, la curva roja de minmod colapsa en ambos extremos y se extiende en un trapecio. El cian THINC-BVD se aferra a la solución exacta gris y conserva su rectángulo. Baja $\beta$ a 1.1 y THINC se ablanda como TVD; cerca de 2.0 es más nítido.

Lo que dejó dudas al reproducirlo#

La promesa es limpia, pero tres cosas trabaron la programación.

La decisión BVD es sensible a la calidad del detector. Confundir un extremo suave (smooth extremum) con una discontinuidad aplana ese pico de forma antinatural. El BVD simple de arriba se protege con el chequeo de monotonía, pero los sistemas multidimensionales necesitan descomposición en variables características.

$\beta$ se mueve con el problema. El valor $\beta = 1.8$ está ajustado a bancos de prueba 1D. En turbulencia compresible fuerte la sobre-compresión puede crear escalones artificiales, así que un $\beta$ adaptativo merece consideración.

La afirmación de costo depende del contexto. "THINC es casi gratis porque las celdas discontinuas son pocas" es cierto si solo se mira la reconstrucción. Pero una prueba BVD con muchas ramas puede disparar divergencia de warps en GPU. Llevarlo a un código FV no estructurado como OpenFOAM implica fundir la detección por celda en el bucle de flujo de caras, lo cual no es tan limpio como una malla estructurada 1D.

Fijar beta en 1.8#

Tras cinco vueltas, THINC-BVD muestra pico 1.00 y ancho de 30 celdas. Minmod, en las mismas condiciones, muestra 0.84 y 23 celdas. Si los extremos suaves quedan prensados en escalones, ese valor fijo sobra. La opción que duda de un solo mando en cada onda es separar THINC de MP5.

Dibujar el salto con tanh, o barrer después un contacto ya difuso, son decisiones distintas. Si la reconstrucción ya sostiene el contacto, un afilado posterior queda ocioso. Si no lo sostiene, conviene mirar el detector antes de subir $\beta$ otra vez.