El mismo resultado con el 18 % de las celdas — criterios de marcado AMR y reflujo coarse-fine

Un búfer más estrecho que el CFL por el intervalo de regrid deja salir la interfaz del parche, y omitir el reflujo cambia el signo de la deriva.

En un cálculo de 50 millones de celdas, ¿cuántas deciden realmente el resultado? Una onda de choque, una capa de cortadura, un frente de llama. Suelen ser un puñado de puntos porcentuales del total. El resto gasta tiempo resolviendo con esmero zonas que ya eran suaves.

Al bajar n_buf a 0 y empujar regrid every hasta 30, crecen celdas rojas a lo largo de la interfaz. El sensor va adelante y la malla no alcanza.

threshold0.35n_buf2max level3regrid every14

leaf 0 / 1 · under-resolved 0

set n_buf to 0 and drag “regrid every” up to 30: red cells appear along the layer every cycle and clear the instant the mesh is rebuilt — that is the feature running out of its own patch. Push n_buf back to 2 and the red stops, but watch the leaf-cell count climb. Raising the threshold thins the patch the other way.

Al subir n_buf a 2 el rojo desaparece y sube el número de celdas hoja. Ese número es el precio de ensanchar el búfer.

El refinamiento adaptativo de malla AMR (adaptive mesh refinement: sembrar celdas solo donde la solución las exige) ataca esa proporción de frente. El sensor que decide dónde subdividir, y la masa que se escapa en la interfaz después de subdividir. Los dos quedan en números.

11 mil celdas de 65 mil#

Conviene empezar contando la ganancia. Se toma un dominio de $256^2$ con una sola capa de cortadura de Kelvin–Helmholtz (la interfaz que se enrolla entre dos capas con diferencia de velocidad). Con malla uniforme son 65 536 celdas. Al envolver esa misma interfaz con AMR de 3 niveles, las celdas hoja quedan en 11 776. Un 18 %.

Esa proporción no es casualidad: la fija la dimensión. En un dominio de $D$ dimensiones, la interfaz tiene $D-1$ . Si la malla uniforme tiene $N_{\rm uni} = (L/h)^D$ celdas, las que cubren solo la interfaz son

N_{\rm AMR} \sim \left(\frac{L}{h}\right)^{D-1} = N_{\rm uni}^{(D-1)/D}

donde $L$ es el tamaño del dominio y $h$ el espaciamiento de la malla más fina. En 2D el exponente vale 1/2; en 3D, 2/3. Es decir: al duplicar la finura, la malla uniforme multiplica sus celdas por 8 y el AMR solo por 4. La ventaja crece con la resolución. Por eso se usa AMR.

Al meter esta malla en la estimación, el número deja de ser vago. La uniforme $256^2$ tiene 65 536 celdas. La suma de hojas que se cuenta más adelante es 11 776, y $11776/65536=0.180$ . Una interfaz recta en 2D sería solo $(L/h)^{1}=256$ celdas, el 0.39 % de la uniforme. La capa enrollada, más el búfer y las hojas gruesas, es lo que llega al 18 %. Si en 3D se parte el espaciamiento a la mitad, la uniforme crece por $2^3=8$ y el AMR de la interfaz por $2^{2}=4$ . El exponente $(D-1)/D$ vale 2/3 en 3D, así que la fracción frente a la malla uniforme sigue bajando al afinar.

El gradiente no sirve para decidir dónde subdividir#

El primer criterio de marcado que se le ocurre a cualquiera es $|\nabla \phi| > \epsilon$ . Falla por una cuestión de unidades. El gradiente de presión va en Pa/m; el de densidad, en kg/m⁴. Hay que recalibrar $\epsilon$ para cada campo, para cada problema e incluso para cada nivel.

El criterio que Löhner publicó en 1987 elimina el problema por normalización: divide la diferencia segunda por la suma de los valores absolutos de las diferencias primeras.

E_i = \sqrt{ \sum_{d} \frac{\left(\phi_{i+e_d} - 2\phi_i + \phi_{i-e_d}\right)^2}{\Big(|\phi_{i+e_d}-\phi_i| + |\phi_i-\phi_{i-e_d}| + \varepsilon\big(|\phi_{i+e_d}| + 2|\phi_i| + |\phi_{i-e_d}|\big)\Big)^2} }

$e_d$ es el vecino en la dirección $d$ y $\varepsilon$ es el filtro de ruido (típicamente 0.01–0.05). Como numerador y denominador tienen la misma dimensión, $E_i$ es adimensional y su valor queda acotado aproximadamente en $[0,1]$ . Un único umbral de 0.3–0.4 funciona igual para la presión que para la densidad, y tanto en el nivel 0 como en el nivel 3.

Un detalle importante: este sensor no detecta rasgos, detecta resolución. Si la interfaz ya está bien resuelta con 4 o 5 celdas en ese nivel, la diferencia segunda se achica y el sensor se apaga. Es una buena propiedad — subdivide solo lo necesario y se detiene solo. Pero si lo que se busca es "todo choque hasta el nivel máximo, sin discusión", este sensor por sí solo no alcanza. En ese caso, en lugar de inflar el término de filtro de ruido $\varepsilon$ , se agrega con un OR un criterio físico aparte (por ejemplo, $\Delta p / p > 0.1$ ).

El búfer es el seguro que cubre hasta el próximo regrid#

Subdividir solo las celdas que encendió el sensor se rompe en el paso siguiente, porque la malla no se rehace en cada paso. El regrid corre una vez cada 4 a 20 pasos, y en ese intervalo la interfaz sigue moviéndose. Por eso las celdas marcadas se engordan en $n_{\rm buf}$ celdas.

El ancho de búfer necesario sale de una cuenta directa. Si en el nivel $\ell$ se convierte a número de celdas la distancia que recorre el rasgo durante $N_{\rm regrid}$ pasos,

n_{\rm buf} \;\ge\; \frac{|u|_{\max}\,\Delta t_\ell\,N_{\rm regrid}}{h_\ell} \;=\; \nu\,N_{\rm regrid}

donde $\nu$ es el número de CFL de ese nivel. Con CFL 0.4 y regrid cada 10 pasos hacen falta al menos 4 celdas. La relación es independiente del nivel: con subciclado, $\Delta t_\ell$ y $h_\ell$ se reducen en la misma proporción.

El piso, escrito, es $n_{\rm buf}\ge 0.4\times 10=4$ . El búfer de 2 del código queda por debajo. Más estrecho que eso, y la interfaz sale de la zona marcada antes del próximo regrid. Esas celdas rojas son lo que muestra el dial del principio con n_buf en 0.

Bloque, celda o parche — lo que define la unidad de corte#

Con el mismo resultado de marcado, la unidad con la que se subdivide cambia tanto el conteo de celdas como la complejidad del código.

Enfoque	Unidad de refinamiento	Estructura de datos	Refinamiento excedente	Implementaciones típicas
Basado en bloques	Bloque de tamaño fijo ( $8^3$ , etc.)	Octree de bloques	Grande	PARAMESH, FLASH
Basado en celdas	Una celda	Árbol por celda	Nulo	OpenFOAM `hexRef8`, RAGE
Basado en parches	Parche rectangular de tamaño libre	Lista de cajas	Pequeño	Chombo, BoxLib/AMReX

El enfoque por bloques tiene la estructura de datos más simple y buena localidad de caché, pero basta una sola celda marcada para que se subdivida el bloque entero. El enfoque por celdas no genera ningún exceso, aunque cada búsqueda de vecinos implica recorrer el árbol. Ahí entra OpenFOAM: dynamicRefineFvMesh es el motor, hexRef8 el cortador (un hexaedro en ocho) y refinementHistory el historial para poder deshacer. Un árbol cuaternario que solo hace crecer las celdas que lo piden usa esa misma unidad. No parte un bloque entero por una sola celda marcada. El enfoque por parches queda en el medio y, dentro del rectángulo, permite conservar los bucles de malla estructurada, lo que favorece la vectorización.

Una misma cara con dos respuestas distintas#

Acá empieza el problema de verdad. Se observa la cara donde se tocan los niveles $\ell$ y $\ell+1$ . Del lado grueso esa cara es una sola; del lado fino son $r^{D-1}$ ( $r$ es la razón de refinamiento, normalmente 2). Y con subciclado, además, los pasos de tiempo difieren.

\Delta t_\ell = \frac{\Delta t_0}{r^{\ell}}

Mientras la malla gruesa avanza un paso, la fina avanza $r$ veces. Por eso el flujo $F^c$ que la malla gruesa calculó a través de esa cara y el flujo que la malla fina realmente empujó, repartido en $r^{D-1}$ caras y $r$ subpasos, son números distintos. La diferencia entre ambos es

\delta F_f = \frac{1}{r^{D}}\sum_{s=1}^{r}\sum_{f'\subset f} F^{s}_{f'} \;-\; F^{c}_{f}

y eso es masa fabricada por la interfaz. Aparece masa que no existía o desaparece la que sí estaba. Ocurre incluso usando un esquema conservativo: la conservación se cumple dentro de una malla, no en el lugar donde dos mallas se encuentran.

Anotar en el libro mayor y saldar de una vez#

La solución es simple. Al arrancar el paso grueso se anota $F^c_f$ en un libro mayor (registro de flujos, flux register). Mientras la malla fina corre sus $r$ subpasos, se acumula en el mismo registro el flujo real. Cuando termina el paso grueso, la diferencia $\delta F_f$ se devuelve a la celda gruesa que está por fuera del parche.

\phi_c \;\leftarrow\; \phi_c \;\mp\; \frac{\Delta t_c}{h_c}\,\delta F_f

El signo lo fija el lado de la cara en que se encuentra esa celda. Las celdas interiores al parche no se tocan: ya fueron sobrescritas con el promedio de los valores finos. Esta única línea devuelve la masa de todo el dominio hasta la precisión de máquina.

Con $r=2$ y $D=2$ , $r^{D}=4$ . Las $r^{D-1}=2$ caras finas y los $r=2$ subciclos se dividen por ese 4 para quedar en la misma escala que un flujo grueso. Sin la división, $\delta F$ sale con el tamaño equivocado. Saldar el registro igual deja la masa descuadrada.

Conteo de celdas y deriva de masa, medidos con código#

Primero el marcado. Se aplica el sensor de Löhner sobre una instantánea de la capa de cortadura, se agrega el búfer, se apilan niveles por bloques y se cuentan las celdas hoja.

import numpy as np
 
N_EFF, BLOCK, MAX_LEVEL = 256, 4, 2   # resolucion equivalente 256^2, 4x4 celdas por bloque, niveles 0~2
EPS_L, THRESH, N_BUF = 0.02, 0.35, 2  # filtro de ruido / umbral de marcado / bufer (en celdas)
 
 
def shear_layer(n, t=1.35):
    """Instantanea del enrollamiento de Kelvin-Helmholtz: sin solver, solo un campo."""
    x = (np.arange(n) + 0.5) / n
    xx, yy = np.meshgrid(x, x, indexing='ij')
    warp = 0.06 * np.sin(2 * np.pi * xx + t) + 0.025 * np.sin(4 * np.pi * xx - 2 * t)
    return np.tanh((yy - 0.5 - warp) / 0.012)
 
 
def shift(a, d, ax):
    """Acceso a vecinos: x periodico, y de gradiente nulo; sin saltos falsos en la pared."""
    return np.roll(a, d, ax) if ax == 0 else np.pad(a, 1, mode='edge')[1:-1, 1 + d:a.shape[1] + 1 + d]
 
 
def lohner_sensor(f):
    """Diferencia segunda normalizada. Al ser adimensional, un solo umbral sirve para todos los niveles."""
    e2 = np.zeros_like(f)
    for ax in (0, 1):
        p, m = shift(f, -1, ax), shift(f, 1, ax)
        num = np.abs(p - 2.0 * f + m)
        den = np.abs(p - f) + np.abs(f - m) + EPS_L * (np.abs(p) + 2 * np.abs(f) + np.abs(m))
        e2 += (num / np.maximum(den, 1e-30)) ** 2
    return np.sqrt(e2)
 
 
def grow(mask, width):
    """Bufer: impide que el rasgo salga del parche antes del proximo regrid."""
    for _ in range(width):
        out = mask.copy()
        for ax in (0, 1):
            out |= shift(mask, 1, ax) | shift(mask, -1, ax)
        mask = out
    return mask
 
 
def tag_blocks(f, thresh, n_buf):
    """Sensor por celda -> bufer por celda -> se refina todo bloque con al menos una celda marcada."""
    tagged = grow(lohner_sensor(f) > thresh, n_buf)
    nb = f.shape[0] // BLOCK
    return tagged.reshape(nb, BLOCK, nb, BLOCK).any(axis=(1, 3))
 
 
def leaf_cells(field):
    """Recorre desde el nivel grueso hacia abajo y cuenta solo las celdas que realmente se resuelven."""
    counts, live = [], None
    for lev in range(MAX_LEVEL + 1):
        n = N_EFF >> (MAX_LEVEL - lev)
        f = field.reshape(n, N_EFF // n, n, N_EFF // n).mean(axis=(1, 3))
        nb = n // BLOCK
        child = np.zeros((nb, nb), bool) if lev == MAX_LEVEL else tag_blocks(f, THRESH, N_BUF)
        live = np.ones((nb, nb), bool) if live is None else live
        counts.append(int((live & ~child).sum()) * BLOCK * BLOCK)
        live = np.kron(live & child, np.ones((2, 2), bool))
    return counts
 
 
counts = leaf_cells(shear_layer(N_EFF))
total, uniform = sum(counts), N_EFF * N_EFF
for lev, c in enumerate(counts):
    print(f'  level {lev}  h = 1/{N_EFF >> (MAX_LEVEL - lev):<3d}  leaf cells = {c:6d}')
print(f'  AMR total     = {total}')
print(f'  uniform 256^2 = {uniform}  ->  {100 * total / uniform:.1f} % of the cells')

La salida es esta.

  level 0  h = 1/64   leaf cells =   3280
  level 1  h = 1/128  leaf cells =   1520
  level 2  h = 1/256  leaf cells =   6976
  AMR total     = 11776
  uniform 256^2 = 65536  ->  18.0 % of the cells

Al cambiar N_BUF a 0, 1, 2 y 4, el total se mueve de 9616 a 10 864, luego 11 776 y finalmente 14 224. El búfer de 2 celdas cuesta 2160 celdas, es decir 3.3 puntos porcentuales respecto de la malla uniforme.

Las 2160 celdas añadidas entre búfer 0 y búfer 2 son cerca del 22 % de 9616. Con 4 celdas el total llega a 14 224, esto es 4608 celdas por encima del caso sin búfer. Por nivel, las 6976 más finas son el 59 % de las 11 776. Las 3280 del nivel 0 cubren la zona suave con un espaciamiento grueso. El 18 % es la suma de esas tres capas, no la interfaz sola.

Ahora el reflujo. Se montan dos niveles sobre una advección escalar 2D, se agrega subciclado y se mide la masa total.

import numpy as np
 
NC, R, CFL, NSTEP = 48, 2, 0.4, 60   # malla gruesa / razon de refinamiento / CFL / pasos gruesos
BOX = (12, 28, 16, 32)               # esquinas del parche (en indices de la malla gruesa)
U, V = 1.0, 0.6
 
 
def upwind_faces(f, dx, dy):
    """Flujo donor-cell en todas las caras x e y del bloque periodico. Devuelve (Fx, Fy)."""
    fx = U * (f if U > 0 else np.roll(f, -1, 0))          # la cara i esta a la izquierda de la celda i
    fy = V * (f if V > 0 else np.roll(f, -1, 1))
    return np.roll(fx, 1, 0), np.roll(fy, 1, 1)
 
 
def march_block(f, fx, fy, dt, dx, dy):
    return f - dt / dx * (np.roll(fx, -1, 0) - fx) - dt / dy * (np.roll(fy, -1, 1) - fy)
 
 
def gaussian_patch(n, x0, y0, s):
    c = (np.arange(n) + 0.5) / n
    xx, yy = np.meshgrid(c, c, indexing='ij')
    return np.exp(-((xx - x0) ** 2 + (yy - y0) ** 2) / s ** 2)
 
 
def two_level_run(reflux):
    i0, i1, j0, j1 = BOX
    dx, dxf = 1.0 / NC, 1.0 / NC / R
    dt = CFL * dx / (abs(U) + abs(V))
    dtf = dt / R
 
    coarse = gaussian_patch(NC, 0.32, 0.42, 0.09)
    fine = np.kron(coarse[i0:i1, j0:j1], np.ones((R, R)))   # el parche arranca consistente con el valor grueso
 
    for _ in range(NSTEP):
        cfx, cfy = upwind_faces(coarse, dx, dx)
        # el flujo que la malla gruesa "creyo" mover a traves del borde del parche
        edge_c = {'lo_x': cfx[i0, j0:j1].copy(), 'hi_x': cfx[i1, j0:j1].copy(),
                  'lo_y': cfy[i0:i1, j0].copy(), 'hi_y': cfy[i0:i1, j1].copy()}
        coarse = march_block(coarse, cfx, cfy, dt, dx, dx)
 
        edge_f = {k: np.zeros_like(v) for k, v in edge_c.items()}
        for _ in range(R):                                   # subciclado: R pasos finos por cada paso grueso
            g = np.zeros((fine.shape[0] + 2, fine.shape[1] + 2))
            g[1:-1, 1:-1] = fine
            g[0, 1:-1] = np.repeat(coarse[i0 - 1, j0:j1], R)  # ghost: valor grueso inyectado como constante a trozos
            g[-1, 1:-1] = np.repeat(coarse[i1, j0:j1], R)
            g[1:-1, 0] = np.repeat(coarse[i0:i1, j0 - 1], R)
            g[1:-1, -1] = np.repeat(coarse[i0:i1, j1], R)
            gfx, gfy = upwind_faces(g, dxf, dxf)
            fine = march_block(g, gfx, gfy, dtf, dxf, dxf)[1:-1, 1:-1]
            for k, s in (('lo_x', gfx[1, 1:-1]), ('hi_x', gfx[-1, 1:-1]),
                         ('lo_y', gfy[1:-1, 1]), ('hi_y', gfy[1:-1, -1])):
                edge_f[k] += s.reshape(-1, R).mean(axis=1) / R   # promedio en la cara y en el tiempo
 
        coarse[i0:i1, j0:j1] = fine.reshape(i1 - i0, R, j1 - j0, R).mean(axis=(1, 3))
 
        if reflux:                                            # liquidacion del libro mayor
            coarse[i0 - 1, j0:j1] -= dt / dx * (edge_f['lo_x'] - edge_c['lo_x'])
            coarse[i1, j0:j1] += dt / dx * (edge_f['hi_x'] - edge_c['hi_x'])
            coarse[i0:i1, j0 - 1] -= dt / dx * (edge_f['lo_y'] - edge_c['lo_y'])
            coarse[i0:i1, j1] += dt / dx * (edge_f['hi_y'] - edge_c['hi_y'])
 
    mask = np.ones((NC, NC), bool)
    mask[i0:i1, j0:j1] = False
    return (coarse * mask).sum() * dx * dx + fine.sum() * dxf * dxf
 
 
m0 = gaussian_patch(NC, 0.32, 0.42, 0.09).sum() / NC ** 2
for tag, on in (('reflux off', False), ('reflux on ', True)):
    m = two_level_run(on)
    print(f'  {tag}:  mass = {m:.12f}   drift = {(m - m0) / m0:+.3e}')

  reflux off:  mass = 0.025702960114   drift = +1.006e-02
  reflux on :  mass = 0.025446894886   drift = +0.000e+00

Al bajar pulse sigma, la desviación roja crece y el $\delta F$ del registro crece en la misma medida.

ratio r4CFL0.60pulse sigma0.050

step 0 · drift 0.00e+0 vs 0.00e+0

the pulse is harmless until it reaches the pink face. From that step on the red trace walks away from zero and never comes back, while the green one stays pinned at machine zero. Sharpen the pulse (sigma down) or raise r and the red excursion grows — the register grows with it, because it is exactly the same number with the opposite sign.

Cuando el pulso toca la cara rosada, solo la traza roja se separa de cero. La verde sigue pegada. La deriva +1.006e-02 sin reflujo es esa línea roja.

Un 1 % en 60 pasos. La resta de las dos masas impresas da cerca de $2.56\times 10^{-4}$ . Dividida por la masa con reflujo, coincide con la deriva +1.006e-02. Y ese valor tampoco crece de forma monótona: a los 15 pasos era −1.3 %, a los 30 pasos −0.77 % y a los 60 pasos +1.0 %. El signo cambia cada vez que el blob entra o sale del parche. Ver una curva así durante un estudio de convergencia lleva a sospechar del esquema, cuando el culpable es una interfaz sin reflujo.

Balanceo de carga — una curva que pone los bloques en fila#

Correr AMR en paralelo trae un problema nuevo: en cada regrid, el número de bloques por procesador cambia. El rango que quedó con la interfaz encima multiplica sus celdas por cinco, mientras que el que solo tiene zona suave se queda igual.

La receta estándar es la curva de llenado de espacio SFC (space-filling curve: una curva que recorre en una sola fila una malla multidimensional). Con una curva de Morton (orden Z) o de Hilbert se asigna un índice unidimensional a cada bloque hoja y esa fila se parte en tantos tramos iguales como rangos haya. Gracias a la localidad de la curva, los índices contiguos suelen ser también vecinos en el espacio, así que una partición equitativa resulta además una partición de poca comunicación. Así trabajan p4est y AMReX. Un particionador de grafos (ParMETIS, Zoltan, Scotch) da mejor calidad de partición, pero hay que volver a ejecutarlo en cada regrid, con un costo muy superior al de la SFC. Como en AMR la malla cambia seguido, la SFC suele ganar.

Sin reflujo, la curva de convergencia señala al esquema#

El marcado se deja adimensional. Una diferencia segunda normalizada permite que un umbral cruce campos y niveles. El gradiente crudo hay que reajustarlo en cada nivel. El búfer se calcula como $\nu N_{\rm regrid}$ . Con CFL 0.4 y 10 pasos son 4 celdas. Poner 1 celda a ojo y hacer regrid cada 20 pasos deja una parte grande del avance con el rasgo ya fuera del parche.

Tratar el reflujo como opcional hace que la cara coarse-fine fabrique masa. El +1 % del paso 60 ni siquiera es monótono. Una deriva que cambia de signo se le atribuye fácil al esquema de advección. El mismo tipo de error aparece cuando correr el refine una línea más tarde infló el error de un ciclo temporal LBM. Si el orden está mal, la respuesta en la interfaz sigue mal aunque las celdas sean pequeñas.