[Reseña] Segundo orden en una sola etapa — un solver GRP para el modelo multifásico de Kapila

Una derivada temporal en la interfaz compra el segundo orden temporal sin otra etapa de Runge-Kutta.

El segundo orden en el tiempo se memorizaba como varias etapas de Runge-Kutta. Este artículo saca ese orden de una sola etapa. La herramienta es un solver GRP (Generalized Riemann Problem) que devuelve el valor de interfaz y su derivada temporal juntos.

Artículo: Tuowei Chen, Zhifang Du, Generalized Riemann Problem Method for the Kapila Model of Compressible Multiphase Flows, arXiv:2310.08241 (2023).

Al apagar la casilla GRP, el control de pendiente se bloquea y el valor de interfaz se congela en el tiempo. Al subir el CFL, solo avanza la marca de medio paso que entra al flujo.

advection speed a: 1.00left slope uₓ: 0.50CFL: 0.80GRP (off = Godunov)

Con la pendiente en cero el punto cian se detiene y coincide con Godunov. Al cambiar el signo de la velocidad de advección $a$ , cambia la celda que aporta el valor. El anillo ámbar es el estado en $\Delta t/2$ .

Segundo orden en una etapa — el problema que resuelve GRP#

Un paso de volumen finito queda fijado por el flujo de interfaz $\mathbf{F}_{i+1/2}$ . Godunov resuelve un problema de Riemann y obtiene un número. Ese número es constante en el tiempo. El segundo orden temporal obliga entonces a repetir el mismo operador espacial, al estilo Runge-Kutta.

GRP evalúa el flujo en el tiempo intermedio $t^{n+1/2}$ .

\mathbf{U}(\mathbf{x}_{i+1/2}, t^{n+1/2}) = \mathbf{U}^{n,*}_{i+1/2} + \frac{\Delta t}{2}\left(\frac{\partial \mathbf{U}}{\partial t}\right)^{n,*}_{i+1/2}

$\mathbf{U}^{n,*}$ es el estado de Riemann en la interfaz. $(\partial \mathbf{U}/\partial t)^{n,*}$ es su derivada temporal instantánea. Sumados como un paso de Taylor, el estado de interfaz avanza $\Delta t/2$ hacia el futuro. Una evaluación.

De constante a tramos a lineal a tramos — Riemann frente a GRP#

Un problema de Riemann tiene datos constantes a tramos a cada lado. Un problema GRP tiene datos lineales a tramos. Las ondas que salen del salto ya no son rectas, y el solver devuelve una derivada temporal junto al estado de Riemann.

En la advección escalar $u_t + a u_x = 0$ , la pendiente de la celda $i$ es $\sigma_i$ . Si $a>0$ , el estado de Riemann en la interfaz es la extrapolación desde la celda izquierda.

U^{*}_{i+1/2} = u_i + \sigma_i \frac{\Delta x}{2}

$u_i$ es el promedio de celda, $\sigma_i$ la pendiente interior, $\Delta x$ el ancho de celda. El punto cian de arriba es esta expresión. Poner la pendiente en cero y caer a Godunov sale de la misma expresión.

El procedimiento de Lax-Wendroff — convertir una derivada temporal en espacial#

La derivada temporal sale de la propia ecuación.

\frac{\partial u}{\partial t} = -a \frac{\partial u}{\partial x} = -a\,\sigma_i

El lado izquierdo pasa a ser la derivada espacial de la derecha. La reconstrucción ya conoce $\sigma_i$ . El valor de interfaz en el medio paso es

U^{n+1/2}_{i+1/2} = U^{*}_{i+1/2} - a\,\sigma_i \frac{\Delta t}{2}

y el flujo es $F = a\,U^{n+1/2}$ . En Euler o Kapila el mismo paso es $\partial u/\partial t = -A\,\partial u/\partial x$ , con jacobiano $A$ . El artículo linealiza las ondas acústicas así y resuelve el GRP de forma aproximada.

Python — midiendo el orden de convergencia de un solver GRP escalar#

Un campo inicial suave $u_0=\sin(2\pi x)$ se advecta un periodo. Se comparan los órdenes $L_1$ de Godunov y GRP en advección lineal periódica.

import numpy as np
 
# Adveccion lineal 1D  u_t + a u_x = 0,  dominio periodico [0,1]
# Godunov de primer orden vs GRP de segundo orden en una etapa (Lax-Wendroff)
 
def cell_gradient(u, dx):
    # pendiente central (sin limitar) para un orden limpio en solucion suave
    return (np.roll(u, -1) - np.roll(u, 1)) / (2.0 * dx)
 
def godunov_step(u, a, dx, dt):
    # constante a tramos: valor de interfaz aguas arriba (a > 0)
    flux = a * u                       # F_{i+1/2} = a u_i
    return u - dt / dx * (flux - np.roll(flux, 1))
 
def grp_step(u, a, dx, dt):
    g = cell_gradient(u, dx)           # pendiente espacial sigma_i
    u_star = u + g * dx / 2.0          # solucion de Riemann en interfaz U*
    u_half = u_star - a * g * dt / 2.0 # medio paso de Lax-Wendroff U^{n+1/2}
    flux = a * u_half                  # F_{i+1/2}^{n+1/2}
    return u - dt / dx * (flux - np.roll(flux, 1))
 
def l1_error(num, exact, dx):
    return np.sum(np.abs(num - exact)) * dx
 
def convergence_study(a=1.0, cfl=0.4, t_end=1.0):
    print(f"{'N':>5} {'Godunov':>16} {'GRP':>16}")
    prev = {}
    for N in [40, 80, 160, 320]:
        dx = 1.0 / N
        x = (np.arange(N) + 0.5) * dx
        u0 = np.sin(2 * np.pi * x)     # condicion inicial suave
        dt = cfl * dx / abs(a)
        nsteps = int(round(t_end / dt))
        dt = t_end / nsteps            # aterrizar exactamente en t_end
        ug = u0.copy(); ur = u0.copy()
        for _ in range(nsteps):
            ug = godunov_step(ug, a, dx, dt)
            ur = grp_step(ur, a, dx, dt)
        exact = np.sin(2 * np.pi * (x - a * t_end))
        eg = l1_error(ug, exact, dx)
        er = l1_error(ur, exact, dx)
        og = f"{np.log2(prev['g']/eg):.2f}" if 'g' in prev else "  - "
        orr = f"{np.log2(prev['r']/er):.2f}" if 'r' in prev else "  - "
        print(f"{N:>5} {eg:.2e}({og}) {er:.2e}({orr})")
        prev = {'g': eg, 'r': er}
 
convergence_study()

La corrida imprime esto.

    N          Godunov              GRP
   40 1.63e-01(  - ) 1.31e-03(  - )
   80 8.76e-02(0.90) 2.69e-04(2.28)
  160 4.54e-02(0.95) 6.32e-05(2.09)
  320 2.31e-02(0.97) 1.55e-05(2.03)

Al pasar el control de $\lambda\Delta t$ de 1.6, el Euler explícito se queda dentro de la banda $[0,1]$ o la rompe. Mover la fracción de volumen inicial y el valor de equilibrio muestra si Crank-Nicolson sigue dentro de esa banda.

stiffness λΔt: 1.60α₀ (initial): 0.85α_eq (target): 0.10

Hasta $\lambda\Delta t=3$ la curva cian baja al equilibrio y la roja oscila. El orden impreso 2.03 pertenece a la advección suave. El control es lo que ocurre cuando esa derivada de medio paso entra en una integral de fuente.

En N=40 el cociente de errores es $1.63\times 10^{-1}/1.31\times 10^{-3}$ , cerca de 124. La frase de cien veces menor nombra esta malla. En N=320 el cociente es $2.31\times 10^{-2}/1.55\times 10^{-5}$ , cerca de 1490. Refinar ocho veces, de 40 a 320, recorta Godunov por 7.1 y GRP por 85. Ocho al cuadrado es 64. Los órdenes 2.28, 2.09 y 2.03 rebasan un poco 64 al asentarse, y 0.90, 0.95, 0.97 suben hacia 1.

El CFL del código es 0.4. La corrección de medio paso mide $|a\,\sigma_i\,\Delta t/2|$ y la extrapolación $|\sigma_i\,\Delta x/2|$ . El cociente es el CFL. El instante del flujo devuelve el 40 por ciento del empuje de la pendiente. Con CFL 0.1 esa parte baja al 10 por ciento. La prueba del orden sigue. Entra menos derivada temporal en un flujo.

En N=80 el cociente ronda 326, y en N=160 ronda 718. Duplicar la malla recorre 124, 326, 718, 1490.

El rígido término fuente de fracción de volumen del modelo de Kapila#

La ecuación de fracción de volumen del modelo de cinco ecuaciones de Kapila, reducción del modelo de siete ecuaciones BN, carga una fuente rígida.

\frac{\partial \alpha_1}{\partial t} + \mathbf{u}\cdot\nabla\alpha_1 = \alpha_1 K\,\nabla\cdot\mathbf{u}

$\alpha_1$ es la fracción de volumen de la fase 1. $K$ es el coeficiente de compactación y expansión tomado de la velocidad del sonido de Wood. Cómo se dobla esa velocidad con la fracción de volumen es la misma pregunta que la velocidad de Wood al caer dentro de una mezcla. La fuente puede empujar $\alpha_1$ fuera de $[0,1]$ . La integración explícita cruza la cota y diverge cuando crece la rigidez.

La derivada temporal de GRP deja conocido el estado de interfaz en el tiempo nuevo $t^{n+1}$ . Una estimación de Gauss-Green da entonces el promedio de celda de $\nabla\cdot\mathbf{u}$ en ese instante. La fuente se integra con Crank-Nicolson, semi-implícito, y cabe escribir un límite de paso que mantenga la fracción de volumen dentro de la cota. Las curvas de $\lambda\Delta t$ de arriba son un modelo de ese límite.

Lo que resultó dudoso al reproducirlo#

El GRP escalar salió de segundo orden. El sistema de Kapila es más áspero. Una pendiente central sin limitar oscila en un choque. El orden apareció porque la prueba era suave. En un salto hace falta un limitador como minmod, y el orden local vuelve a uno.

La velocidad del sonido de Wood no es monótona en la fracción de volumen. Por eso el número de Mach se mueve dentro del emborronado numérico de una interfaz. La derivada temporal de GRP sigue la compactación y la expansión. No quita la falta de monotonía. El artículo rodea eso con robustez, no con un arreglo directo.

Los códigos de producción suelen juntar la división de Strang con Runge-Kutta, así que un GRP de una etapa no entra tal cual. Conservar la ventaja exige reescribir el flujo en forma Lax-Wendroff. El porte sigue costando. El paso acústico linealizado se remonta al GRP de dinámica de gases de Ben-Artzi y Falcovitz.

Esperar el orden con la fuente todavía fuera del flujo#

El 2.03 impreso es un número para una pendiente suave. Con un limitador encendido en un choque, ese lugar vuelve localmente al primer orden. Antes de borrar etapas porque la tabla dice 2, tiene que existir el presupuesto para reescribir el flujo.

Dejar la fuente fuera del flujo hace que otro integrador temporal herede la misma rigidez. La relajación que lleva dos presiones y luego las funde en una o bien guarda la cota de la fracción de volumen en ese lado separado, o bien vuelve a cortar $\Delta t$ . El 2.03 impreso no responde eso.