[Reseña] Elegir la reconstrucción con menor salto

Un salto puro de una celda anula una diferencia, así que THINC queda fuera antes de la comparación de saltos.

En 2018 el grupo de Xiao en Tokyo Tech midió un espesor que pocos quieren publicar. Flujo bifásico compresible, 1000 pasos. El salto que empezó con una celda se había extendido a ocho. WENO tampoco lo detuvo. La disipación de un solo paso puede parecer pequeña. El tiempo igual arruina la interfaz.

La reconstrucción que trajeron es MUSCL-THINC-BVD. Cada celda carga dos candidatas. Se queda la de menor salto en las caras.

Al mover el deslizador de β de 1.4 a 2.2, la tanh dentro de una celda de salto se estrecha, y los hombros de la gaussiana responden con dientes.

N = 80β = 1.60CFL = 0.40

MUSCL THINC BVDt = 0.00THINC picked: 0% of cells

Conviene mirar si las barras cian se encienden solo en los bordes de la onda cuadrada. También si esas barras dejan de multiplicarse cerca de 1.6.

Resumen en una página#

Autores / revista: Deng, Inaba, Xie, Shyue, Xiao. Journal of Computational Physics 371 (2018) 945–966.
Problema: En el modelo de cinco ecuaciones para flujo bifásico compresible, las interfaces materiales se extienden con el tiempo.
Propuesta: Construir MUSCL (tramos suaves) y THINC (saltos) en cada celda, y conservar el salto de cara más pequeño. Esa regla es BVD, boundary variation diminishing.
Qué cambia: Sin antidifusión de posproceso y sin compresión artificial. La misma regla cae sobre la fracción de volumen y el resto de variables conservadas, así que la coherencia entre ellas viene sola.

Dos exigencias en conflicto#

La numérica bifásica compresible quiere dos cosas a la vez. Las zonas suaves piden exactitud y poca disipación. Los saltos, o sea interfaces y choques, piden monotonía y un perfil delgado. Una sola función rara vez hace las dos.

MUSCL es una recta dentro de cada celda. La monotonía se sostiene. La pendiente llega solo a primer orden. Cada cruce de la interfaz recorta un poco. En 1000 pasos, una celda pasa a ocho. THINC ajusta una tanh dentro de la celda y encierra el escalón en una o dos celdas. Poner ese ajuste sobre una zona suave fabrica escalones falsos.

El artículo arranca con una decisión. No mezclar las dos en una función. Elegir una por celda.

MUSCL — el difusor confiable#

La candidata de base es MUSCL con limitador de pendiente minmod.

\tilde{q}_i^{\,\text{MUSCL}}(x) = \bar{q}_i + \sigma_i (x - x_i), \quad \sigma_i = \text{minmod}(\bar{q}_i - \bar{q}_{i-1},\; \bar{q}_{i+1} - \bar{q}_i)

$\bar q_i$ es el promedio de celda. $\sigma_i$ es la pendiente en el centro. minmod devuelve la magnitud menor cuando las dos diferencias comparten signo, y cero cuando no.

Los valores de cara se leen así.

q^{L,\text{MUSCL}}_{i+1/2} = \bar{q}_i + \tfrac{1}{2}\sigma_i, \quad q^{R,\text{MUSCL}}_{i-1/2} = \bar{q}_i - \tfrac{1}{2}\sigma_i

No hay oscilaciones. En un salto, solo crece el espesor, y lo hace despacio.

THINC — una función monótona que imita saltos#

THINC, tangente de hipérbola para capturar interfaces, ajusta una tangente hiperbólica dentro de la celda.

\tilde{q}_i^{\,\text{THINC}}(x) = \bar{q}_{\min} + \tfrac{\bar{q}_{\max}}{2}\left[1 + \theta \tanh\!\left(\beta\!\left(\tfrac{x - x_{i-1/2}}{\Delta x} - \tilde{x}_i\right)\right)\right]

$\bar q_{\min}$ y $\bar q_{\max}$ son el mínimo y la amplitud tomados de las vecinas. $\theta$ es la dirección del salto. $\beta$ es el espesor. $\tilde x_i$ es el centro del salto, resuelto para que el promedio de celda no se mueva.

$\beta$ entre 1.4 y 2.0 es la banda estable. 1.6 es el valor de costumbre. Valores mayores encierran el salto en una celda, lo cual se ve limpio, y cortan dientes en las zonas suaves. Un registro que fija la tanh del lado del choque, en lugar de elegir candidata, está en la reconstrucción que apila TENO sobre THINC. Allí las dos se mezclan por pesos.

BVD — medir la variación entre dos candidatas#

Con dos candidatas, se mide cuál salto de cara es menor cuando se encuentran con las reconstrucciones vecinas. La variación total de frontera de la celda $i$ se define así.

\text{TBV}_i^{P} = |q^{L,\text{MUSCL}}_{i-1/2} - q^{R,P}_{i-1/2}| + |q^{L,P}_{i+1/2} - q^{R,\text{MUSCL}}_{i+1/2}|

$P$ es MUSCL o THINC. Las dos vecinas se quedan en MUSCL. Solo cambia la celda central.

La regla de selección es esta.

\tilde{q}_i^{\,\text{BVD}}(x) = \begin{cases} \tilde{q}_i^{\,\text{THINC}}(x) & \text{if } \delta < C_i < 1 - \delta,\ (\bar{q}_{i+1}-\bar{q}_i)(\bar{q}_i-\bar{q}_{i-1}) > 0,\ \text{TBV}^{\text{THINC}}_i < \text{TBV}^{\text{MUSCL}}_i \\ \tilde{q}_i^{\,\text{MUSCL}}(x) & \text{otherwise} \end{cases}

$C_i$ dice dónde se sienta el salto dentro de la celda. $\delta$ vale cerca de $10^{-4}$ .

La primera compuerta multiplica las dos diferencias de un lado y exige un producto estrictamente positivo. Un factor nulo la anula. Un escalón puro de 0 a 1, de una celda, tiene esa forma. La celda anterior al salto tiene diferencia a la derecha y ninguna a la izquierda. Se vuelve a MUSCL antes de armar la variación de frontera.

La segunda compuerta es la posición. $\delta$ queda cerca de $10^{-4}$ . Un promedio dentro de esa fracción de la amplitud vecina no es elegible para THINC. Las crestas planas salen.

Consistencia sobre el modelo de cinco ecuaciones#

El artículo no se detiene en un concurso de reconstrucciones. El modelo de cinco ecuaciones lleva la fracción de volumen $\alpha_1$ , las densidades de fase $\alpha_1\rho_1$ y $\alpha_2\rho_2$ , el momento $\rho u$ y la energía total $E$ . La misma regla BVD se aplica a las cinco.

Los métodos anteriores afilaban solo la fracción de volumen y luego parcheaban el resto para que la presión no sonara en la interfaz. BVD obliga a las cinco variables a compartir una decisión en cada celda. Una celda de salto usa THINC en las cinco. Una celda suave usa MUSCL en las cinco. La coherencia viene sola. No hace falta una etapa aparte de antidifusión.

El árbol de decisión BVD en NumPy#

import numpy as np
 
def muscl_minmod_edges(q):
    """ Recibe (q[i-1], q[i], q[i+1]) y devuelve los valores izquierdo y derecho de MUSCL """
    qm, q0, qp = q
    a, b = q0 - qm, qp - q0
    if a * b <= 0.0:
        slope = 0.0
    else:
        slope = np.sign(a) * min(abs(a), abs(b))
    return q0 - 0.5 * slope, q0 + 0.5 * slope  # qL, qR
 
def thinc_jump_edges(q, beta=1.6, eps=1e-20):
    """ Salto ajustado con tanh: valores de borde de THINC """
    qm, q0, qp = q
    qmin = min(qm, qp)
    qmax = max(qm, qp) - qmin
    if qmax < 1e-12:
        return q0, q0
    theta = np.sign(qp - qm)
    C = (q0 - qmin + eps) / (qmax + eps)
    if C <= 1e-6 or C >= 1 - 1e-6:
        return q0, q0
    B = np.exp(theta * beta * (2.0 * C - 1.0))
    A = (B / np.cosh(beta) - 1.0) / np.tanh(beta)
    qR = qmin + 0.5 * qmax * (1.0 + theta * A)
    num = 1.0 + theta * A * np.tanh(beta) + theta * np.tanh(beta)
    den = 1.0 + A * np.tanh(beta)
    qL = qmin + 0.5 * qmax * (num / den)
    return qL, qR
 
def bvd_decide(q_window, beta=1.6, delta=1e-4):
    """
    q_window: arreglo de longitud 5 [q[i-2], q[i-1], q[i], q[i+1], q[i+2]]
    devuelve: (qL_i, qR_i, picked_thinc)
    """
    qm2, qm1, q0, qp1, qp2 = q_window
    monotone = (qp1 - q0) * (q0 - qm1) > 0.0
    qmin = min(qm1, qp1)
    qmax = max(qm1, qp1) - qmin
    C = 0.5 if qmax < 1e-12 else (q0 - qmin) / qmax
    eligible = monotone and (delta < C < 1.0 - delta)
 
    qL_M, qR_M = muscl_minmod_edges([qm1, q0, qp1])
    if not eligible:
        return qL_M, qR_M, False
 
    qL_T, qR_T = thinc_jump_edges([qm1, q0, qp1], beta=beta)
    # Comparación de TBV con vecinas fijadas en MUSCL
    _, qR_left = muscl_minmod_edges([qm2, qm1, q0])
    qL_right, _ = muscl_minmod_edges([q0, qp1, qp2])
    tbv_M = abs(qR_left - qL_M) + abs(qR_M - qL_right)
    tbv_T = abs(qR_left - qL_T) + abs(qR_T - qL_right)
    if tbv_T < tbv_M:
        return qL_T, qR_T, True
    return qL_M, qR_M, False
 
# Prueba: perfil suave + salto
N = 80
x = (np.arange(N) + 0.5) / N
q = np.where((x >= 0.6) & (x <= 0.8), 1.0, 0.0)
q += 0.5 * np.exp(-((x - 0.25) / 0.04) ** 2)
 
picks = 0
for i in range(2, N - 2):
    _, _, used_thinc = bvd_decide(q[i-2:i+3])
    picks += int(used_thinc)
print(f"Tasa de elección de THINC: {picks}/{N-4}")
# Esperado: solo unas pocas celdas alrededor del salto

El comentario del fragmento cita 4 de 76 celdas con THINC. La onda del código es 0 o 1, así que una diferencia es cero y el producto falla. Esas cuatro celdas piden una rampa del mismo signo en dos o tres celdas. Tras 1000 pasos, ocho celdas ya son esa rampa. Con una celda de espesor, esta rama no llama a THINC.

En una recta las diferencias vecinas coinciden, la pendiente minmod las iguala y la variación de frontera de MUSCL es cero. La tanh de THINC abre una holgura en la cara y la desigualdad se queda con MUSCL. Una gaussiana más ancha que la celda es esa recta en lo local. Por eso los tramos suaves rechazan el cambio.

Una gaussiana suave está localmente cerca de una recta, así que el cambio se rechaza. Una onda cuadrada que sigue teniendo una celda de espesor tampoco abre la compuerta.

Quién gana en cada celda — verlo en directo#

La advección de arriba empuja tres reconstrucciones a la vez y pinta una barra en cada celda donde BVD elige THINC. Subir $\beta$ de 1.4 a 2.2 estrecha solo esas celdas. Demasiado grande, y aparecen dientes en los hombros de la gaussiana. El número de barras deja de crecer cerca de 1.6, y por eso ese valor se quedó. El cian aparece en los bordes solo después de que la onda cuadrada se ha emborronado un poco. Mientras el espesor es de una celda, las barras siguen vacías.

Hasta dónde se le puede creer#

El orden observado se queda por debajo de dos. Una candidata es MUSCL, así que la exactitud suave se detiene donde se detiene MUSCL. Cambiar a una candidata de orden alto, como WENO, puede llegar a quinto orden. Este artículo eligió la versión simple.

TBV supone que las vecinas son MUSCL. Un caso en el que ambos lados quieren THINC a la vez se ignora dentro de un solo barrido. La sección 4 del artículo sigue viéndose aceptable bajo ese atajo. Una forma más general quedó para después.

En varias dimensiones, $\beta$ se parte según la normal de la cara. La normal de la interfaz se estima con el algoritmo de Young. En esquinas y quiebres la normal falla, y el salto se inclina hacia un lado.

La integración en el tiempo va aparte. El artículo usa SSPRK3. La decisión BVD tiene que rehacerse en cada etapa de Runge–Kutta. Si la decisión parpadea entre etapas, la interfaz tiembla. Un trabajo posterior volvió a ese temblor.

Donde la elección se detiene en un salto de una celda#

Mientras la interfaz sigue teniendo una celda de espesor, esta rama nunca llama a THINC. El producto de diferencias se vuelve positivo solo después de que el espesor se ha extendido sobre dos o tres celdas. Retener una interfaz de una celda y poder elegir solo después de que ya se extendió son dos tareas que se desencuentran.

Si ese desencuentro se filtra al límite del paso de tiempo es la pregunta siguiente. Una interfaz delgada sigue teniendo un techo de advección. La advección de interfaz que se quedó en CFL 0.05 después de retirar el término capilar pregunta dónde se sienta ese techo. Si BVD de verdad encierra el salto en una celda, ¿hasta dónde puede crecer el paso que mueve esa celda?

[Reseña] Elegir la reconstrucción con menor salto — Deng (2018) MUSCL-THINC-BVD