Advectar sin dividir direcciones — El esquema upwind no dividido de Bell-Colella-Glaz

Tras una vuelta el error de masa es redondeo en ambos esquemas, pero los picos quedan en 0.349 y 0.806. Aprobar conservación no significa que el disco siga ahí.

A finales de los años 80, John Bell y Phillip Colella, en Lawrence Berkeley, notaron una cicatriz extraña en sus simulaciones de explosiones. Un frente de llama colocado en diagonal respecto a los ejes de la malla se iba difuminando y torciendo en esa misma diagonal con el tiempo. El culpable no era la física, sino la discretización. La división direccional (dimensional splitting) —resolver una dirección a la vez— distorsionaba la advección diagonal. La salida es el esquema upwind no dividido de segundo orden de Bell-Colella-Glaz (BCG). Al dividir direcciones, un disco en rotación se vuelve un rombo. Un término de derivada transversal tapa ese hueco.

Al subir el CFL de 0 a 0.9, el punto ámbar se acerca al promedio de celda. Con el limiter en none, la recta junto al salto pasa los valores vecinos.

CFL = 0.50

Cyan = limited PLM slope inside each cell · Amber dot = predicted left edge state at n+1/2 (advecting velocity u>0). Pick none near the jump to see the slope overshoot the neighbor values.

En 0.9 el peso de la pendiente casi desaparece. Al volver a minmod, ese exceso se recorta.

La ecuación es la advección conservativa más simple que existe. Un campo de velocidad sin divergencia (divergence-free) $\mathbf{u}$ transporta un escalar $s$ .

\frac{\partial s}{\partial t} + \nabla\cdot(\mathbf{u}\,s) = 0

Aquí $s$ es un escalar pasivo como concentración o temperatura, y $\mathbf{u}$ es un campo de velocidad dado. Al actualizar los promedios de celda con volúmenes finitos (FV) se obtiene esta forma.

s_{ij}^{n+1} = s_{ij}^{n} - \frac{\Delta t}{\Delta x}\left(F^x_{i+\frac12,j} - F^x_{i-\frac12,j}\right) - \frac{\Delta t}{\Delta y}\left(F^y_{i,j+\frac12} - F^y_{i,j-\frac12}\right)

$F^x = u\,s^{n+1/2}$ es un flujo centrado en el tiempo (time-centered) en la cara. Toda la precisión se reduce a qué tan bien se estima el valor en la cara $s^{n+1/2}$ .

Dividir una dirección a la vez tuerce la rotación#

La división direccional parte un problema 2D en una actualización 1D en x seguida de otra 1D en y. La implementación es fácil: basta un único solver 1D. El problema es el flujo diagonal.

Cuando un escalar se mueve a $45^\circ$ , resolver x primero y luego y da una respuesta distinta que el orden inverso. Ambos pierden el término cruzado (cross term) que vive entre los dos pasos 1D. En un flujo en rotación ese error se acumula en cada paso, y un disco redondo se deforma en un rombo. Alternar el orden con división de Strang solo cancela la asimetría direccional del error de segundo orden; la causa raíz permanece.

La receta de BCG es sencilla: resolver x e y a la vez. Dentro de un mismo paso, predecir todos los valores de cara al mismo instante $t^{n+1/2}$ e incorporar de antemano el cambio transversal a esa predicción.

Convertir tiempo en espacio dentro del predictor#

La idea central es empujar el valor de cara medio paso adelante en tiempo y en espacio con una expansión de Taylor. En la cara derecha de la celda $i$ , el estado izquierdo se lee así.

s^{n+1/2}_{i+\frac12,\,L} = s^n_{ij} + \frac{\Delta x}{2}\,\partial_x s + \frac{\Delta t}{2}\,\partial_t s

Dejar la derivada temporal tal cual no sirve de nada. Se sustituye la ecuación de advección $\partial_t s = -u\,\partial_x s$ para cambiar la derivada temporal por una espacial (el truco de Cauchy-Kowalewski).

s^{n+1/2}_{i+\frac12,\,L} = s^n_{ij} + \frac{1}{2}\left(1 - \frac{\Delta t}{\Delta x}\,u_{i+\frac12,j}\right)\Delta_x s_{ij}

$\Delta_x s_{ij}$ es la pendiente limitada (limited slope) de la celda $i$ , y $u_{i+1/2,j}$ es la velocidad advectiva en la cara. El factor $(1 - \mathrm{CFL})$ entre paréntesis carga el avance temporal. A medida que CFL se acerca a 1, la contribución del estado izquierdo se reduce; en 0 queda una simple extrapolación espacial.

El estado derecho, visto desde la celda vecina, se construye de forma simétrica.

s^{n+1/2}_{i+\frac12,\,R} = s^n_{i+1,j} - \frac{1}{2}\left(1 + \frac{\Delta t}{\Delta x}\,u_{i+\frac12,j}\right)\Delta_x s_{i+1,j}

Elegir entre los dos candidatos es la regla upwind de Godunov. El signo de la velocidad de cara decide el viento.

s^{n+1/2}_{i+\frac12} = \begin{cases} s_{i+\frac12,\,L}, & u_{i+\frac12} > 0 \\[2pt] s_{i+\frac12,\,R}, & u_{i+\frac12} < 0 \end{cases}

En la figura de arriba, un CFL de 0.9 acerca el punto ámbar al promedio, y none deja que la recta pase los vecinos. Esa diferencia es el factor de avance temporal del predictor de justo arriba.

La derivada transversal cose las esquinas#

BCG inyecta el acoplamiento que la división deja caer directamente en el predictor. Al estado izquierdo de una cara en x se le resta medio paso de la divergencia del flujo transversal.

s^{n+1/2}_{i+\frac12,\,L} \;\mathrel{-}=\; \frac{\Delta t}{2}\,\frac{\partial (v\,s)}{\partial y}\bigg|_{ij}

$v$ es la velocidad en y, y $\partial_y(v s)$ es el cambio advectivo transversal en la celda $ij$ . Para estimar este término primero hace falta obtener los estados predichos 1D $\hat s^y$ sobre las caras en y. Así que BCG tiene dos pasadas. Primero predecir el estado upwind 1D en cada dirección, luego reutilizar esa predicción en la corrección transversal, y por último elegir el valor de cara final con la misma regla upwind. Este "acoplamiento de esquina (corner coupling)" transporta la advección diagonal sin torsión.

Ponerle riendas a la pendiente — minmod y superbee#

La precisión del predictor sale de la pendiente $\Delta_x s$ . Usar la pendiente centrada en crudo es de segundo orden, pero sobrepasa (overshoot) cerca de las discontinuidades. Un limitador de pendiente (slope limiter) le pone riendas. El más conservador, minmod, toma la menor de las dos diferencias laterales, o cero cuando difieren en signo.

\Delta_x s_{ij} = \operatorname{minmod}\!\left(s_{ij}-s_{i-1,j},\; s_{i+1,j}-s_{ij}\right)

Superbee mantiene la misma monotonía pero exprime la pendiente lo más empinada posible, afilando más los frentes de contacto. En la demo 1D de arriba, al alternar entre minmod y superbee se observa que la recta PLM del lado superbee es más empinada.

Poner la misma rienda sobre el flujo, en una dimensión, es el limitador TVD sobre Lax–Wendroff. BCG cuelga el limitador de la pendiente de celda, no del flujo.

Girando un disco en numpy#

Ahora corremos el BCG 2D no dividido. El problema de juguete es rotación de cuerpo rígido (solid-body rotation). Se flota un disco escalar en un campo sin divergencia que gira alrededor del centro, y se compara el upwind de primer orden contra BCG en la misma malla.

import numpy as np
 
N = 64
dx = 1.0 / N
omega = 2 * np.pi                      # una vuelta por unidad de tiempo
 
def face_velocities():
    j = np.arange(N)
    i = np.arange(N)
    u = -omega * ((j + 0.5) * dx - 0.5)        # velocidad en cara x (depende de y)
    u = np.broadcast_to(u[None, :], (N, N)).copy()
    v = omega * ((i + 0.5) * dx - 0.5)         # velocidad en cara y (depende de x)
    v = np.broadcast_to(v[:, None], (N, N)).copy()
    return u, v
 
def minmod_pair(a, b):                          # pendiente que preserva monotonia
    return np.where(a * b > 0, np.where(np.abs(a) < np.abs(b), a, b), 0.0)
 
def limited_slopes(s):
    slx = minmod_pair(s - np.roll(s, 1, 0), np.roll(s, -1, 0) - s)
    sly = minmod_pair(s - np.roll(s, 1, 1), np.roll(s, -1, 1) - s)
    return slx, sly
 
def donor_fluxes(s, u, v):                       # upwind de primer orden (celda donante)
    Fx = np.where(u > 0, u * np.roll(s, 1, 0), u * s)
    Fy = np.where(v > 0, v * np.roll(s, 1, 1), v * s)
    return Fx, Fy
 
def bcg_fluxes(s, u, v, dt):                     # predictor no dividido de segundo orden
    slx, sly = limited_slopes(s)
    # 1) estados predichos 1D en la direccion normal
    sLx = np.roll(s, 1, 0) + 0.5 * (1 - dt/dx * u) * np.roll(slx, 1, 0)
    sRx = s - 0.5 * (1 + dt/dx * u) * slx
    hatX = np.where(u > 0, sLx, np.where(u < 0, sRx, 0.5 * (sLx + sRx)))
    sBy = np.roll(s, 1, 1) + 0.5 * (1 - dt/dx * v) * np.roll(sly, 1, 1)
    sTy = s - 0.5 * (1 + dt/dx * v) * sly
    hatY = np.where(v > 0, sBy, np.where(v < 0, sTy, 0.5 * (sBy + sTy)))
    # 2) correccion transversal (de esquina) — reusar los estados predichos
    divVy = (np.roll(v * hatY, -1, 1) - v * hatY) / dx
    divUx = (np.roll(u * hatX, -1, 0) - u * hatX) / dx
    sLx -= 0.5 * dt * np.roll(divVy, 1, 0)
    sRx -= 0.5 * dt * divVy
    sBy -= 0.5 * dt * np.roll(divUx, 1, 1)
    sTy -= 0.5 * dt * divUx
    # 3) elegir el valor de cara final con la misma regla upwind
    sX = np.where(u > 0, sLx, np.where(u < 0, sRx, 0.5 * (sLx + sRx)))
    sY = np.where(v > 0, sBy, np.where(v < 0, sTy, 0.5 * (sBy + sTy)))
    return u * sX, v * sY
 
def advance(s, u, v, dt, scheme):
    Fx, Fy = bcg_fluxes(s, u, v, dt) if scheme == 'bcg' else donor_fluxes(s, u, v)
    return s - dt * ((np.roll(Fx, -1, 0) - Fx) / dx + (np.roll(Fy, -1, 1) - Fy) / dx)
 
def make_disk():
    X, Y = np.meshgrid((np.arange(N)+0.5)*dx, (np.arange(N)+0.5)*dx, indexing='ij')
    return ((X - 0.5)**2 + (Y - 0.78)**2 < 0.13**2).astype(float)
 
u, v = face_velocities()
Umax = omega * 0.5 * np.sqrt(2)
dt = 0.6 * dx / Umax                    # CFL = 0.6
steps = int(round(1.0 / dt))            # exactamente una vuelta
 
for scheme in ('upwind1', 'bcg'):
    s = make_disk()
    mass0 = s.sum()
    for _ in range(steps):
        s = advance(s, u, v, dt, scheme)
    print(f"{scheme:8s}  peak={s.max():.3f}  "
          f"mass_err={abs(s.sum()-mass0)/mass0:.2e}")

La salida se ve así.

upwind1   peak=0.349  mass_err=3.1e-16
bcg       peak=0.806  mass_err=2.8e-16

Al pasar el esquema de 1st-order upwind a BCG unsplit, el pico recorre el camino de 0.349 a 0.806, los mismos números de la salida.

CFL = 0.60

A scalar disk rotates rigidly about the center. Switch to 1st-order upwind and watch the disk melt into a ring within one revolution; BCG unsplit keeps its edge far longer.

En 1st-order upwind el disco se abre en un anillo antes de completar una vuelta, con el borde borroso. En BCG unsplit la frontera dura mucho más. Con CFL 0.9, BCG también se vuelve algo más rugoso. La brecha frente al primer orden se mantiene.

Ambos esquemas conservan la masa hasta la precisión de máquina. Lo garantiza la diferenciación conservativa de flujos. Lo que cambia es el pico. El upwind de primer orden aplasta el disco a menos de la mitad de su altura en una vuelta, mientras que BCG conserva el 80% de la altura original.

Qué queda después de dos vueltas#

El código ya fija una vuelta. N es 64, así que el espaciado es 1/64. La velocidad de esquina es 2π por 0.5 por la raíz de 2, cerca de 4.44. Con CFL 0.6 el paso es cerca de 0.00211, y el tiempo 1 lleva 474 pasos. El centro del disco (0.5, 0.78) queda a 0.28 de (0.5, 0.5) y se mueve a cerca de 1.76. El paso usa la velocidad de esquina 4.44, no 1.76.

El diámetro 0.26 son 16.6 celdas, unas 217 celdas de escalar. Con CFL 0.6 el peso izquierdo de la pendiente es 0.2 y el derecho 0.8. En el tope 0.8 el peso izquierdo baja a 0.1. Los picos 0.349 y 0.806 son esos pesos más el término transversal, apilados 474 veces.

Los errores de masa 3.1e-16 y 2.8e-16 son redondeo. La altura perdida es 0.651 contra 0.194, cerca de 3.4 veces más en el upwind que en BCG. El término transversal resta cerca de 0.00106 por la divergencia transversal, una vez por paso, 474 veces por vuelta. El difuminado lo decide cuánta pendiente queda en una celda, no el número 64.

Dónde se filtra al portarlo#

El CFL es más estricto que en un esquema dividido. El 2D no dividido propaga información por la diagonal, así que su límite de estabilidad es menor que el 1D. En la práctica conviene empezar con $\mathrm{CFL}\lesssim 0.8$ .
Guardar siempre la velocidad en las caras. Interpolar la velocidad centrada en celda a las caras rompe la propiedad de divergencia nula, y la masa que se creía conservada se fuga. La disposición MAC (escalonada) es la respuesta.
No descartar el término transversal. Si solo se implementa la predicción normal, queda un MUSCL con división dimensional sin más. El código corre, pero la torsión rotacional permanece.
El limiter se anula en los extremos. En la cima de un montículo la pendiente se recorta a cero y cae a primer orden localmente. Por eso el pico se hunde poco a poco. Para preservarlo más, conviene revisar limiters que conservan extremos.

La elección en la que la masa cuadra y el disco igual muere#

Un pico de 0.349 es lo que el upwind de primer orden deja tras una vuelta. BCG queda en 0.806. Los errores de masa son 3.1e-16 y 2.8e-16, así que un registro de conservación aprueba los dos. Al mirar el disco, solo uno sigue siendo un disco.

Si el upwind de primer orden sirve para un escalar en diagonal, se decide por esa brecha. Cuando la altura del contacto se puede tirar, la celda donante alcanza. Conservar la altura exige meter el término transversal en el mismo paso. Al salir de una malla cartesiana, minmod no viaja gratis: el trato de los extremos cambia. El limitador de Barth–Jespersen muestra cómo un extremo falso se corta con el mínimo y el máximo de los vecinos. Esa es la decisión que sigue.