ρ y u coinciden hasta el último dígito, pero la deformación se infla un 45 % — el reescalado de no equilibrio en el refinamiento de malla LBM

Lo único que cruza intacto la frontera de nivel son ρ y u. La parte de no equilibrio hay que reescribirla multiplicada por τ_f/(m·τ_c).

ρ y u coincidían hasta el último dígito, y aun así el esfuerzo se infló#

Al acoplar refinamiento adaptativo de malla (AMR) a un código de Boltzmann en red (LBM) aparece un objeto nuevo: la frontera de niveles. Ahí se superponen una capa de malla gruesa y una de malla fina. En esa capa hay que trasladar los valores de un lado al otro.

Para validar ese traslado suele mirarse la densidad y la velocidad. Y con eso solo, el error pasa desapercibido. Aunque las funciones de distribución $f_i$ se copien tal cual al nivel fino, $\rho$ y $\mathbf{u}$ siguen coincidiendo hasta el último decimal. Lo que no coincide es el esfuerzo. Este texto precisa por cuánto se desvía y por qué justamente por ese factor.

La conclusión primero: el factor es el inverso de $\tilde\tau_f/(m\tilde\tau_c)$ . En una configuración habitual da 1.45; en configuraciones de alto Reynolds se acerca a 2.

Lo que debe seguir igual al cruzar de nivel y lo que está obligado a cambiar#

El principio de la conversión entre niveles es uno solo: las cantidades conservadas deben coincidir. $\rho$ , $\mathbf{u}$ y $p$ son magnitudes físicas y no tienen motivo para cambiar porque cambie la malla.

Sea $m = \Delta x_c/\Delta x_f$ la razón de refinamiento. El subíndice $c$ indica malla gruesa y $f$ , malla fina. Con escalado acústico (acoustic scaling), $\Delta t$ se reduce en la misma proporción.

\Delta x_f = \frac{\Delta x_c}{m}, \qquad \Delta t_f = \frac{\Delta t_c}{m}

Así, la velocidad en unidades de red $\hat{u} = u\,\Delta t/\Delta x$ es la misma en ambos niveles. Como $\rho$ también lo es, la distribución de equilibrio $f_i^{\mathrm{eq}}(\rho, \hat{u})$ resulta exactamente el mismo valor en los dos niveles. Lo único que hay que modificar es la desviación respecto al equilibrio, es decir $f_i^{\mathrm{neq}} = f_i - f_i^{\mathrm{eq}}$ .

Conviene manipularlo directamente en la simulación siguiente.

tau_c0.80m2

Both lattices start from the same sine and advance over the same physical time. With tau rescaled the orange curve sits on the blue one; press tau copied and the fine level suddenly relaxes a different fluid. At step 0 the amplitudes are 1.0000 and 1.0000; the measured flux ratio is 0.0000 against the predicted r = 0.6875. Drag tau_c towards 0.51 and r falls to 1/2 — that is where copying a population costs you a factor of two.

Dos mallas reales resuelven en paralelo la misma región física durante el mismo tiempo físico. Al mover el deslizador tau_c, vale la pena vigilar si measured q_f/q_c en la escalera de la derecha se mantiene pegado a rescale r. Con el botón tau copied, el nivel fino pasa a ser un fluido de otra viscosidad y la curva naranja se despega de la azul.

Si se fija la viscosidad, τ se mueve detrás#

La viscosidad física no puede cambiar porque cambie el nivel. Escrita con el tiempo de relajación en unidades de red $\tilde\tau$ , queda así.

\nu = c_s^2\left(\tilde\tau - \tfrac{1}{2}\right)\frac{\Delta x^2}{\Delta t}

$c_s^2$ es el cuadrado de la velocidad del sonido de red y $\Delta x^2/\Delta t$ tiene unidades de difusividad. Ese $-\tfrac{1}{2}$ es el rastro que dejó la integración trapezoidal, tratado aparte en el Δt/2 que deja la discretización LBM.

Bajo escalado acústico, $\Delta x^2/\Delta t$ queda multiplicado por $1/m$ . Al imponer $\nu_f = \nu_c$ , $\tilde\tau$ sale despejado.

\tilde\tau_f = \frac{1}{2} + m\left(\tilde\tau_c - \frac{1}{2}\right)

Con $m=2$ y $\tilde\tau_c = 0.8$ resulta $\tilde\tau_f = 1.1$ . Lo importante es que no es el doble. Multiplicar $\tilde\tau$ por $m$ sin más descuadra $\nu$ , porque $\tfrac{1}{2}$ no es objeto del escalado.

A la parte de no equilibrio no se le pega solo τ, también Δt#

Si el término de primer orden del desarrollo de Chapman–Enskog se escribe con la aproximación de Grad, la parte de no equilibrio queda proporcional a la tasa de deformación.

f_i^{\mathrm{neq}} = -\frac{w_i \rho \tilde\tau}{c_s^2}\left(\mathbf{c}_i\mathbf{c}_i - c_s^2\mathbf{I}\right) : \hat{\mathbf{S}}

$w_i$ son los pesos, $\mathbf{c}_i$ las velocidades de red y $\hat{\mathbf{S}}$ la tasa de deformación en unidades de red. Su relación con la tasa física $\mathbf{S}$ es $\hat{\mathbf{S}} = \mathbf{S}\,\Delta t$ . Como la deformación física en un mismo punto no depende del nivel, $f^{\mathrm{neq}}$ resulta proporcional a $\tilde\tau\,\Delta t$ .

r \equiv \frac{f_i^{\mathrm{neq},f}}{f_i^{\mathrm{neq},c}} = \frac{\tilde\tau_f\,\Delta t_f}{\tilde\tau_c\,\Delta t_c} = \frac{\tilde\tau_f}{m\,\tilde\tau_c}

Ese $r$ es el factor de reescalado, y coincide con la forma que ordenaron Dupuis y Chopard. Con $\tilde\tau_c = 0.8$ y $m=2$ se obtiene $r = 1.1/1.6 = 0.6875$ : el no equilibrio del nivel fino es el 69 % del que tiene el nivel grueso.

Tres formas de traspaso puestas en la misma tabla#

En la práctica, pasar valores a través de la frontera de niveles se reparte en tres caminos.

Forma	Qué se traspasa	Información necesaria	Costo	Dónde falla
Copia completa	$f_i$ tal cual	Ninguna	El más bajo	$\rho,\mathbf{u}$ salen bien y el esfuerzo se infla $1/r$ veces
Macroscópicas + reconstrucción de Grad	Interpolar $\rho, \mathbf{u}$ y recalcular $f^{\mathrm{neq}}$	Gradiente de velocidad	Intermedio	Hay que volver a estimar el gradiente por diferencias finitas
Interpolación del no equilibrio + escala	$f^{\mathrm{eq}}$ recalculado, $f^{\mathrm{neq}}$ por $r$	$\tilde\tau_c, \tilde\tau_f, m$	Casi el más bajo	Es fácil olvidar la $m$ dentro de $r$

La segunda y la tercera deberían dar el mismo resultado. La tercera es el estándar práctico porque evita recalcular gradientes: con dos $\tilde\tau$ y una $m$ ya sale el factor.

La interpolación espacial en sí no tiene misterio. De gruesa a fina, interpolación trilineal en tres dimensiones; de fina a gruesa, el promedio de $2^d$ celdas, donde $d$ es el número de dimensiones. Lo difícil no es interpolar, sino decidir a qué se le multiplica $r$ después de interpolar.

La deformación recuperada, medida en Python#

Sobre una red D2Q9 se fija una tasa de deformación física y se construyen los $f^{\mathrm{neq}}$ de cada nivel. Luego se recupera la deformación con el $\tilde\tau_f$ del nivel fino y se ponen lado a lado el valor reescalado y el simplemente copiado.

CS2 = 1.0 / 3.0
EX = [0, 1, 0, -1, 0, 1, -1, -1, 1]
EY = [0, 0, 1, 0, -1, 1, 1, -1, -1]
W = [4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36]
 
 
def equilibrium(rho, ux, uy):
    """Equilibrio D2Q9. Con el mismo rho y u, da el mismo valor aunque cambie el nivel."""
    out = []
    u2 = ux*ux + uy*uy
    for i in range(9):
        cu = EX[i]*ux + EY[i]*uy
        out.append(rho*W[i]*(1 + cu/CS2 + cu*cu/(2*CS2*CS2) - u2/(2*CS2)))
    return out
 
 
def grad_neq(rho, tau, sxy):
    """Parte de no equilibrio según Grad. sxy es la deformación cortante en unidades de red del nivel."""
    out = []
    for i in range(9):
        q_xy = EX[i]*EY[i]            # componente xy de Q_i (la diagonal no se empareja con sxy)
        out.append(-(W[i]*rho*tau/CS2) * 2.0 * q_xy * sxy)
    return out
 
 
def recover_strain(fneq, rho, tau, dt):
    """Recupera la deformación cortante en unidades físicas desde el momento de no equilibrio."""
    pi_xy = sum(fneq[i]*EX[i]*EY[i] for i in range(9))
    return -pi_xy / (2.0*rho*CS2*tau*dt)
 
 
def tau_on_level(tau_c, m):
    """Tiempo de relajación que debe tener el nivel fino al fijar la viscosidad."""
    return 0.5 + m*(tau_c - 0.5)
 
 
def nu_physical(tau, dx, dt):
    return CS2*(tau - 0.5)*dx*dx/dt
 
 
rho, ux, uy = 1.0, 0.05, 0.0
s_phys = 0.004          # deformación cortante física — valor independiente del nivel
m = 2                   # razón de refinamiento
dx_c, dt_c = 1.0, 1.0
dx_f, dt_f = dx_c/m, dt_c/m
 
feq_c = equilibrium(rho, ux, uy)
feq_f = equilibrium(rho, ux, uy)
print("max |feq_c - feq_f| = %.3e" % max(abs(a-b) for a, b in zip(feq_c, feq_f)))
print()
print("%5s  %6s  %8s  %8s  %7s  %11s  %14s  %6s" % (
    "tau_c", "tau_f", "nu_c", "nu_f", "r_meas", "tf/(m*tc)", "copied/true", "err%"))
for tau_c in [0.51, 0.55, 0.60, 0.80, 1.20, 2.00]:
    tau_f = tau_on_level(tau_c, m)
    fneq_c = grad_neq(rho, tau_c, s_phys*dt_c)
    fneq_f = grad_neq(rho, tau_f, s_phys*dt_f)
    r = fneq_f[5]/fneq_c[5]
    s_ok = recover_strain(fneq_f, rho, tau_f, dt_f)
    s_bad = recover_strain(fneq_c, rho, tau_f, dt_f)
    print("%5.2f  %6.3f  %8.5f  %8.5f  %7.4f  %11.4f  %14.4f  %6.1f" % (
        tau_c, tau_f, nu_physical(tau_c, dx_c, dt_c), nu_physical(tau_f, dx_f, dt_f),
        r, tau_f/(m*tau_c), s_bad/s_ok, (s_bad/s_ok - 1)*100))
 
print()
worst = grad_neq(rho, 0.51, s_phys)
print("mass moment of f^neq     = %.3e" % sum(worst))
print("momentum moments of f^neq = %.3e, %.3e" % (sum(worst[i]*EX[i] for i in range(9)),
                                            sum(worst[i]*EY[i] for i in range(9))))

max |feq_c - feq_f| = 0.000e+00
 
tau_c   tau_f      nu_c      nu_f   r_meas    tf/(m*tc)     copied/true    err%
 0.51   0.520   0.00333   0.00333   0.5098       0.5098          1.9615    96.2
 0.55   0.600   0.01667   0.01667   0.5455       0.5455          1.8333    83.3
 0.60   0.700   0.03333   0.03333   0.5833       0.5833          1.7143    71.4
 0.80   1.100   0.10000   0.10000   0.6875       0.6875          1.4545    45.5
 1.20   1.900   0.23333   0.23333   0.7917       0.7917          1.2632    26.3
 2.00   3.500   0.50000   0.50000   0.8750       0.8750          1.1429    14.3
 
mass moment of f^neq     = 0.000e+00
momentum moments of f^neq = 0.000e+00, 0.000e+00

Salen tres cosas de una vez. La distribución de equilibrio es idéntica en los dos niveles. Las columnas $\nu_c$ y $\nu_f$ coinciden hasta el último dígito. Y el $r$ medido reproduce la forma cerrada $\tilde\tau_f/(m\tilde\tau_c)$ .

Las dos últimas líneas son el origen del título. Los momentos de orden cero y uno de $f^{\mathrm{neq}}$ son exactamente nulos. Se copie o se reescale, $\rho$ y la cantidad de movimiento no se mueven. Lo único que sale mal es el momento de segundo orden.

Cuanto más se pega τ a 0.5, más se duplica el precio de copiar#

Leer la columna de error de arriba hacia abajo revela la tendencia. Con $\tilde\tau_c = 2.0$ es 14 %. Con $\tilde\tau_c = 0.51$ es 96 %. En forma cerrada aparece el rango de $r$ .

r = \frac{1/2 + m(\tilde\tau_c - 1/2)}{m\,\tilde\tau_c}

Si $\tilde\tau_c \to 1/2$ entonces $r \to 1/2$ ; si $\tilde\tau_c \to \infty$ entonces $r \to 1$ . Con viscosidad grande, copiar casi no deja rastro. Con viscosidad pequeña, el esfuerzo se duplica.

El problema es que el motivo para poner AMR suele ser justamente el alto Reynolds, y ahí $\tilde\tau$ se baja hasta rozar 0.5. La región donde el error de copia estalla con más fuerza es la misma región donde se quiere usar AMR. Los síntomas también confunden: como la masa y la cantidad de movimiento se conservan, nada diverge y solo queda una capa fina de vorticidad a lo largo de la línea de nivel.

Subir $m$ lo empeora. Con $m=4$ y $\tilde\tau_c = 0.8$ resulta $\tilde\tau_f = 1.7$ , de modo que $r = 1.7/3.2 = 0.531$ y el error de copia salta al 88 %.

Agregar un nivel se cobra como m^(d+1)#

El reescalado ocurre solo dos veces: en la explosion, que baja del nivel grueso al fino, y en la coalescence, que sube del fino al grueso. Entre ambas todo es collide-and-stream corriente.

El reloj de abajo permite recorrer un ciclo paso a paso.

m2dwell42

Use next phase to walk the cycle one gate at a time: explosion, 2 fine sub-steps, coalescence. The two dashed lines are the only moments populations cross levels — everything between them is ordinary collide-and-stream. Raise m or switch to d = 3 and the work factor climbs as m^(d+1) = 8; currently at cycle 0, phase explosion.

Con next phase se confirma el orden explosion → $m$ subpasos → coalescence. Las líneas punteadas verde y morada marcan los dos únicos instantes en que las distribuciones cruzan de nivel. Al subir m y d, el work factor del registro de la derecha crece como $m^{d+1}$ .

Ese exponente es la restricción real del diseño AMR. En tres dimensiones con $m=2$ , un parche cuesta 16 veces más. Por eso dónde y cuánto nivel fino se coloca pesa más en el rendimiento que la elección del esquema.

Entonces, ¿por qué la capa de solape debe ser de una sola celda?#

La capa de solape es la capa de celdas donde existen a la vez los nodos de ambos niveles. ¿Por qué una sola? Porque el streaming avanza una casilla por paso. Cuando el nivel fino da un subpaso, las distribuciones que entran desde afuera vienen exactamente de una casilla más allá. Basta con rellenar esa casilla.

Este enfoque es el mismo de las distribuciones que pierde un nodo de frontera. Sea una pared o una frontera de niveles, lo primero es contar "qué $f_i$ quedan vacías tras el streaming". En la pared la respuesta la fija la geometría; en la frontera de niveles, la razón de refinamiento $m$ .

Con una malla centrada en celdas (cell-centered), la capa de solape se maneja mejor que con una centrada en nodos (node-centered). Al no haber nodos superpuestos la propiedad queda clara, y en la partición paralela el destinatario de la comunicación se reduce a una sola lista de celdas. A cambio, en la interpolación de gruesa a fina las posiciones se desfasan media casilla y el esténcil debe ajustarse a ello.

Ya sea escribiendo una interfaz AMR desde cero o leyendo el código ajeno, el diagnóstico cabe en tres líneas. ¿Es $\tilde\tau_f$ igual a $\tfrac{1}{2} + m(\tilde\tau_c - \tfrac{1}{2})$ ? ¿Se multiplica $f^{\mathrm{neq}}$ por $r$ ? ¿Está la $m$ en el denominador de $r$ ? La tercera es la que falta más seguido.