AMR con Quadtree — Dejar que la malla crezca solo donde hace falta

Al pasar tau de 0.05 a 0.4 las colas dejan de partirse; sin la regla 2:1, el flujo de la cara es lo primero que falla.

En 1984, en Stanford, Marsha Berger y Joseph Oliger perseguían una onda de choque. La malla tenía cerca de un millón de celdas. Solo alrededor del 1%, junto al choque, hacía trabajo útil. El 99% restante recalculaba un flujo libre uniforme con la misma precisión. El artículo de ese año en Journal of Computational Physics es el AMR de Berger–Oliger. El mismo hábito sigue dentro de dynamicRefineFvMesh de OpenFOAM, AMReX, Chombo y p4est.

Al arrastrar el punto blanco, la malla se parte alrededor de ese punto; el control threshold cambia la pendiente escalada $\tau$ que decide el corte.

gradient threshold: 0.180max level: 5

drag the white dot to move the feature. higher threshold = coarser grid. cell count vs uniform refinement is the AMR savings.

Al bajar $\tau$ a 0.05, las hojas siguen la cola de la gaussiana. Por encima de 0.4 queda una malla gruesa casi uniforme. Con add feature aparece un segundo pico y solo esos dos núcleos se profundizan; la lectura de la esquina superior izquierda es el conteo frente a una malla uniforme.

Por qué una malla uniforme desperdicia el 90%#

Quien haya corrido un RANS compresible 3D conoce el intercambio. Duplicar el espaciado deja un paso 8 veces más barato, porque $2^3=8$ . Choques, capas de cortadura y separación se emborronan con esa decisión. Refinar el dominio entero para esa banda delgada encoge también las celdas del flujo libre. Memoria y tiempo estallan juntos.

Basta con refinar la banda. El cociente es esta estimación.

\frac{N_\text{AMR}}{N_\text{uniform}} \approx \frac{V_\text{feature}}{V_\text{domain}} \cdot 4^{L} + (1 - \tfrac{V_\text{feature}}{V_\text{domain}})

$L$ es el nivel máximo de refinamiento. $V_\text{feature}$ es el volumen que pide espaciado fino. $V_\text{domain}$ es el dominio entero. En 2D, $4^{L}$ es lo que hace un nivel al conteo: cuatro hijos. Con $L=5$ , $4^5=1024$ . Si el 1% del arranque se toma como fracción de volumen, el conteo frente a la malla más gruesa es $0.01\times 1024+0.99=11.23$ veces mayor. Una malla uniforme de ese mismo nivel es 1024 veces la gruesa. Al dividir ambas, el AMR queda cerca del 1.1% de la malla fina uniforme. Un choque con $V_\text{feature}/V_\text{domain}\sim h$ empuja esa fracción todavía más abajo. El porcentaje de un solo dígito sale de esa división. En 3D, partir una arista por la mitad multiplica las celdas por 8. Repartir ese 8 por todo el dominio carga la factura entera. Dejarlo dentro de la banda detiene el factor 8 ahí.

h-, p-, r- — Tres sabores de adaptación#

Una malla adaptativa — la que cambia con la solución — se separa según qué se mueve.

Adaptación h: reducir el tamaño $h$ de la celda. Es el sentido habitual de AMR.
Adaptación p: dejar las celdas y subir el orden polinómico $p$ . Habitual en DGM y SEM (spectral element method).
Adaptación r (relocalización): fijar el número de celdas y el orden, y deslizar los nodos hacia las zonas de gradiente alto. La memoria queda plana. La distorsión no.

El resto de esta nota se queda en h-AMR, y en el quadtree 2D. En 3D el mismo árbol es un octree.

Quadtree — celdas que poseen sus hijos#

Una malla uniforme cabe en un arreglo. Una malla adaptativa necesita un árbol. Un refine crea cuatro hijos, y cada hijo apunta a su padre. Si un hijo se parte, aparecen nietos.

class QuadCell:
    __slots__ = ("cx", "cy", "size", "level", "children", "parent", "data")
 
    def __init__(self, cx, cy, size, level=0, parent=None):
        self.cx = cx           # x del centro de la celda
        self.cy = cy           # y del centro de la celda
        self.size = size       # longitud del lado
        self.level = level     # raíz = 0, refinar suma 1
        self.children = None   # None si es hoja, lista de 4 si no
        self.parent = parent
        self.data = 0.0        # cantidad conservada (densidad, etc.)
 
    def is_leaf(self):
        return self.children is None
 
    def split(self):
        h = self.size * 0.5
        q = self.size * 0.25
        offsets = ((-q, -q), (+q, -q), (-q, +q), (+q, +q))
        self.children = [
            QuadCell(self.cx + dx, self.cy + dy, h, self.level + 1, self)
            for dx, dy in offsets
        ]
        for ch in self.children:
            ch.data = self.data    # inicializar con el padre (conservativo)

Copiar el valor del padre en split es por conservación. La integral del padre $U \cdot V$ tiene que igualar la suma de los cuatro hijos. Cada hijo ocupa un cuarto del volumen del padre. Con la copia directa, las cuatro piezas suman $4 \times (U \cdot V/4) = U \cdot V$ . Esa copia es solo interpolación de primer orden. Los códigos de producción suelen añadir la pendiente del padre, al estilo MUSCL, y tratar el resultado como prolongación de segundo orden. La pendiente puede entrar. La suma de los cuatro sigue teniendo que coincidir con $U \cdot V$ . Si no coincide, la partición crea o destruye masa.

Dónde cortar — el criterio de refinamiento#

Un quadtree es una caja vacía. El criterio decide qué celda se corta. Cinco formas aparecen todo el tiempo.

Criterio	Forma	Lo que captura
Primera derivada	$h \lvert \nabla u \rvert$	Choques, capas de cortadura
Segunda derivada	$h^2 \lvert \nabla^2 u \rvert$	Regiones de alta curvatura
Frobenius del Hessiano	$\sqrt{\sum H_{ij}^2}\, h^2$	Detección de anisotropía
Estimador de Richardson	$\lvert u_h - u_{2h} \rvert$	Error real de discretización
Umbral físico	$\rho > \rho_, \lvert \omega \rvert > \omega_$	Núcleos de vórtices

dynamicRefineFvMesh de OpenFOAM aplica un umbral al campo que eligió quien usa el código. AMReX deja en blanco un booleano por celda, tagcells(), para que lo complete quien llama. Lo que se etiqueta cambia el árbol más que el contenedor. Seguir con la misma respuesta mientras baja el conteo de celdas continúa en etiquetado y reflujo coarse-fine.

La malla del navegador usa solo la primera fila. Diferencia central de $|\nabla u|$ en el centro de la celda, por el tamaño $h$ . Por encima de $\tau$ , la celda se parte. Si no, queda como hoja.

Equilibrio 2:1 y nodos colgantes#

La cara donde un hijo se encuentra con un vecino es la incómoda. Un lado es una celda grande. El otro lado son dos celdas pequeñas. Esa cara cría un nodo colgante (hanging node). En volumen finito, un nodo colgante rompe el balance de flujo.

Hay dos salidas.

Forzar equilibrio 2:1. Las hojas vecinas pueden diferir como mucho en un nivel. p4est, AMReX y los códigos de la familia Berger hacen esto.
Interpolación non-conforming. Construir el valor de la cara gruesa a partir de las dos celdas pequeñas. Habitual en códigos tipo DGM.

El bucle que impone la regla es corto. Recorrer los candidatos a refinar. Si una hoja vecina es dos o más niveles más gruesa, entra en la cola. Repetir hasta vaciarla.

def enforce_2to1(roots, refine_queue):
    """asegurar que cada celda en refine_queue cumple 2:1 con sus vecinas"""
    pending = list(refine_queue)
    queued = set(id(c) for c in pending)
    while pending:
        cell = pending.pop()
        for nb in face_neighbors(cell, roots):
            if nb.is_leaf() and cell.level - nb.level >= 1:
                if id(nb) not in queued:
                    pending.append(nb)
                    queued.add(id(nb))
    return [c for c in (refine_queue + pending) if c.is_leaf()]

face_neighbors devuelve las hojas que comparten una cara. El recorrido es $O(\log N)$ .

Una malla adaptativa 1D en Python#

Un listado completo del quadtree 2D se hace largo. La misma idea cabe en un árbol binario 1D para un paso. Se planta un pulso gaussiano. Solo se bisecan las celdas que superan el umbral.

import numpy as np
 
def gaussian_pulse(x, x0=0.6, sigma=0.04):
    return np.exp(-((x - x0) ** 2) / (2 * sigma ** 2))
 
def refine_indicator(cell_lo, cell_hi, fn, h_scale=1.0):
    """primera derivada escalada = diferencia absoluta a través de la celda"""
    u_lo = fn(cell_lo)
    u_hi = fn(cell_hi)
    return abs(u_hi - u_lo) * h_scale
 
def adaptive_bisect_1d(fn, tau=0.05, max_level=8):
    """bisecar [0, 1] de forma recursiva según el umbral"""
    leaves = []
    def recurse(lo, hi, level):
        ind = refine_indicator(lo, hi, fn)
        if level >= max_level or ind < tau:
            leaves.append((lo, hi, level))
            return
        mid = 0.5 * (lo + hi)
        recurse(lo, mid, level + 1)
        recurse(mid, hi, level + 1)
    recurse(0.0, 1.0, 0)
    return leaves
 
cells = adaptive_bisect_1d(gaussian_pulse, tau=0.05, max_level=7)
print(f"leaf cells: {len(cells)}")
print(f"vs uniforme nivel 7: {len(cells) / 2**7 * 100:.1f}%")
# salida ejemplo: leaf cells: 41, vs uniforme nivel 7: 32.0%

Una gaussiana ocupa una franja delgada, así que la corrida se detiene en 41 hojas: 32.0% de una línea uniforme de nivel 7. Con dos pulsos, solo los dos núcleos se profundizan. El hueco queda grueso.

Entre 41 hojas y 128 celdas#

Dividir el 41 impreso entre $2^7=128$ da 0.3203. El 32.0% del comentario es ese número, cortado a un decimal en porcentaje. Al restar quedan 87 celdas que nunca se construyen con el espaciado más fino, cerca del 68% de la línea uniforme de nivel 7. En el código, $\sigma=0.04$ es el 4% del intervalo unidad. Solo un entorno de ese pulso supera $\tau=0.05$ . El flujo libre plano casi no tiene diferencia entre los dos extremos, así que la bisección se detiene pronto. Un corte 1D produce dos hijos. Un corte 2D produce cuatro, de modo que la misma franja delgada cae más rápido cuando entra $4^{L}$ de la estimación. Subir $\tau$ de 0.05 a 0.4 hace que las colas que estaban entre esos dos números dejen de partirse.

Trampa del paralelo — load balancing y SFC#

La memoria distribuida choca con otro muro. Un rango que se queda con todas las ramas profundas fija el wall clock de todos. El conteo de celdas dejó de estar equilibrado.

Hay que alinear el árbol en 1D y cortar esa línea. Una curva que llena el espacio (space-filling curve, SFC) — Morton (Z-order) o Hilbert — mantiene cerca en la línea a las celdas que ya son vecinas en el espacio. El índice Morton de p4est intercala los bits de las coordenadas.

def morton_2d(ix, iy, level):
    """intercalar coordenadas enteras (ix, iy) en un código Morton"""
    code = 0
    for b in range(level):
        code |= ((ix >> b) & 1) << (2 * b)
        code |= ((iy >> b) & 1) << (2 * b + 1)
    return code
 
# Z-order sobre una malla 4x4
for iy in range(4):
    for ix in range(4):
        print(morton_2d(ix, iy, 2), end=" ")
    print()
# 0  1  4  5
# 2  3  6  7
# 8  9 12 13
# 10 11 14 15

Ordenar por índice Morton y luego partir la lista entre rangos junta a los vecinos y encoge los mensajes. Los dominios cuadrados siguen saltando en las esquinas. Una curva de Hilbert evita el salto y cuesta más de indexar, así que Morton es el estándar de facto en la práctica.

Marcar el campo equivocado deja caro incluso ese 32%#

La forma del árbol está casi cerrada. Cuatro hijos, y un corte Morton para que los vecinos compartan rango. En un bloque 4×4 el índice empieza en 0, 1, 4, 5 porque los bits de las coordenadas van intercalados. El gasto aparece después. Marcar la densidad mientras el error vive en una capa de cortadura apila esas 41 hojas de nivel 7 en la banda equivocada. El 32% pertenece a un criterio apuntado a $x_0=0.6$ con $\sigma=0.04$ . Cambiar el campo etiquetado mueve $V_\text{feature}$ en la estimación, incluido el factor 11.23 frente a la malla gruesa.

Saltarse el 2:1 es una decisión real. Limitar la diferencia de nivel a 1 deja al nodo colgante con dos celdas pequeñas como máximo en una cara. Quitar el límite obliga a escribir interpolación non-conforming en cada cara. El flujo de la cara se separa entonces de la integral del padre $U\cdot V$ . El árbol puede estar cerrado y la pregunta de tiempo sigue abierta: en qué línea del avance temporal cae el refine. Cuánto mueve el error una línea mal puesta es el tema de el ciclo temporal recursivo en LBM AMR.