[Reseña de paper] Terminar la simulación antes de la tormenta — Casulli–Walters (2000), aguas someras semi-implícitas en mallas ortogonales no estructuradas

La celeridad sale del límite de estabilidad solo mientras el excedente diagonal de la matriz queda en uno.

En junio de 1990 el laboratorio BAW de Hamburgo tenía que reproducir 13 días de marea en el Jade–Weser. 150 000 triángulos. Paso de 5 minutos. Más rápido que el tiempo real, o el siguiente pronóstico no arranca.

Un paso explícito deja que la celeridad de superficie libre $c=\sqrt{gh}$ fije el CFL. A 28 m, $c \approx 16{,}6$ m/s. Con el lado más corto cerca de $\sqrt{10}$ m, el paso permitido es 0,2 s. Trece días son 1 123 200 s. A 0,2 s salen 5 616 000 pasos, los 5,6 millones del artículo. A 5 minutos salen 3 744 pasos. El cociente es exactamente 1 500, porque 300 / 0,2 no depende de qué quincena se elija.

Casulli y Walters (2000) sostienen ese cociente dejando solo el gradiente de superficie libre en el tiempo nuevo. La advección se queda en el tiempo viejo.

Al subir el deslizador de CFL por encima de 1 la curva explícita roja estalla, y al subir θ solo baja la amplitud cian.

CFL ratio (Δt / Δt_explicit): 2.00θ (implicitness): 0.60

CFL > 1 ⇒ explicit blows up; semi-implicit stays bounded for any CFL when θ ≥ 1/2.

En CFL=6 la cresta cian sigue ahí. Al llevar θ hasta 1 esa cresta muere rápido; la franja donde la amplitud aguanta en la práctica va de 0,55 a 0,6.

Metadatos del paper#

Autores: Vincenzo Casulli (Univ. de Trento), Roy A. Walters (USGS)
Revista: International Journal for Numerical Methods in Fluids, 32, 331–348 (2000)
Título: "An unstructured grid, three-dimensional model based on the shallow water equations"
Palabras clave: semi-implicit, shallow water, unstructured orthogonal grid, wetting/drying

Por qué el solver explícito falla en estuarios#

A 10 m, $c \approx 10$ m/s. A 100 m, $c \approx 31$ m/s. La velocidad $u$ suele quedar por debajo de 1 m/s. La marea cruza un tramo en $L/u$ , pero el paso explícito queda atado a $\Delta x/c$ . Con $c/u \approx 10\sim30$ , el código corta el tiempo bastante más fino de lo que el agua pide.

Una malla no estructurada carga todo el límite sobre el lado más corto. Big Lost River mide 1,2 km × 1 km, 12 622 triángulos, lado mínimo cerca de 5 m, límite explícito de unos 0,5 s. Dividir 5 m entre 0,5 s da 10 m/s. No es la celeridad de 16 m/s del canal profundo la que detiene el reloj. Es el lado corto en agua somera.

Los términos lentos siguen explícitos. El rápido pasa a implícito. Cuando solo el gradiente de presión de superficie libre es implícito, $\sqrt{gh}$ sale de la cota de estabilidad. Lo que queda es un límite desde $u$ .

Mallas ortogonales no estructuradas — el dúo Voronoi–Delaunay#

Las diferencias finitas, puestas tal cual sobre una malla no estructurada, se rompen en dos sitios. El espaciado cambia con la forma y el gradiente se inclina. Escapar por coordenadas curvilíneas añade términos de transformación que recortan precisión y estabilidad a la vez.

Casulli y Walters restringieron la malla. Malla ortogonal no estructurada: el segmento entre centros vecinos es perpendicular a la arista compartida. Un polígono de Voronoi y su triangulación de Delaunay dual cumplen eso cuando todo triángulo es agudo. Un ángulo obtuso y el circuncentro sale del triángulo. La unión de centros puede seguir perpendicular a la arista y aun así no cruzar el interior de la celda. Cuánto tiempo la inserción mantiene vacío el circuncírculo es cuánto tiempo sigue válida esta diferencia.

La ortogonalidad compra una diferencia. Con solo la velocidad normal $u_j$ , $\partial \eta / \partial n$ es $(\eta_{i_2} - \eta_{i_1})/\delta_j$ . Sin corrección de curvatura.

División semi-implícita: el método θ separa dos escalas temporales#

El método θ parte el tiempo con un factor de implicitud $\theta \in [1/2, 1]$ . El momento ve el gradiente de superficie libre como $\theta \eta^{n+1} + (1-\theta)\eta^n$ . La ecuación de superficie libre ve la velocidad como $\theta u^{n+1} + (1-\theta)u^n$ .

Momento discretizado, con viscosidad vertical y esfuerzo del viento:

\mathbf{A}_j^n \mathbf{U}_j^{n+1} = \mathbf{G}_j^n - \theta g \frac{\Delta t}{\delta_j}\!\left[\eta_{i(j,2)}^{n+1} - \eta_{i(j,1)}^{n+1}\right]\Delta\mathbf{Z}_j^n

$\mathbf{U}_j^{n+1}$ es la velocidad normal en la arista $j$ . $\mathbf{A}_j^n$ es la tridiagonal de viscosidad vertical, viento y fricción de fondo. $\mathbf{G}_j^n$ guarda advección, Coriolis y fricción horizontal. Cada arista es un sistema $N_z \times N_z$ , y hay $N_s$ de ellos.

La ecuación de superficie libre:

P_i \eta_i^{n+1} = P_i \eta_i^n - \theta \Delta t \sum_{l=1}^{S_i} s_{i,l}\lambda_{j(i,l)} [\Delta\mathbf{Z}^\top \mathbf{U}]_{j(i,l)}^{n+1} - (1-\theta)\Delta t \sum_{l=1}^{S_i} s_{i,l}\lambda_{j(i,l)} [\Delta\mathbf{Z}^\top \mathbf{U}]_{j(i,l)}^{n}

$P_i$ es el área de la celda, $\lambda_j$ la longitud de la arista, $s_{i,l}$ el signo hacia afuera.

Reducción a ecuación de onda y la recompensa SPD#

Acopladas, las dos ecuaciones tienen $N_z N_s + N_p$ incógnitas. Se despeja $\mathbf{U}_j^{n+1}$ en el momento como función de $\eta^{n+1}$ y se sustituye. Si $\mathbf{A}$ es definida positiva, también lo es su inversa.

Lo que cae es una ecuación de onda discreta.

P_i \eta_i^{n+1} - g \Delta t^2 \theta^2 \sum_{l=1}^{S_i} \frac{s_{i,l}\lambda_{j(i,l)}}{\delta_{j(i,l)}}\!\left[(\Delta\mathbf{Z})^\top \mathbf{A}^{-1}\Delta\mathbf{Z}\right]_{j(i,l)}^{n}\!\!\left(\eta^{n+1}_{i[j(i,l),2]} - \eta^{n+1}_{i[j(i,l),1]}\right) = R_i

Una incógnita, $\eta^{n+1}$ . Tamaño $N_p$ . Simétrica, con diagonal fuertemente dominante, definida positiva. PCG encaja. El artículo reporta convergencia en decenas de iteraciones incluso en mallas grandes.

En el juguete 1D no hay estructura vertical, así que $\mathbf{A}$ vale 1 y el $k = g H \Delta t^2 \theta^2 / \Delta x^2$ del código es el coeficiente del vecino. A cuatro veces el límite explícito, $c^2 = gH$ fuerza $k = 16\theta^2$ . Con θ=0,6, $k=5{,}76$ . La diagonal interior $1+2k$ vale 12,52. Las dos entradas fuera de la diagonal suman 11,52. El excedente es siempre 1, y ese 1 es el término de almacenamiento $P_i \eta$ . Al crecer el paso, $k$ se pega a la diagonal y a la fuera de diagonal por igual, y el excedente no crece. En el CFL=6 del deslizador y θ=0,6, $k=12{,}96$ , la diagonal es 26,92 y la suma fuera de diagonal es 25,92. Excedente: sigue en 1. Un paso más implícito no vuelve la matriz más segura. Lo que muere es la amplitud.

Con $\eta$ conocido, $\mathbf{U}_j^{n+1}$ vuelve a una tridiagonal pequeña en cada arista. La velocidad vertical se recupera de la continuidad discreta.

Una versión 1D en Python#

La misma partición. $N=200$ , $L=1$ , $H=1$ , $g=9{,}81$ , θ=0,6.

import numpy as np
from scipy.linalg import solve_banded
 
class CasulliShallowWater1D:
    """Aguas someras 1D lineales — método θ semi-implícito."""
 
    def __init__(self, N=200, L=1.0, H=1.0, g=9.81, theta=0.6):
        self.N, self.L, self.H, self.g, self.theta = N, L, H, g, theta
        self.dx = L / N
        self.eta = np.zeros(N)            # superficie libre en centros de celda
        self.u = np.zeros(N + 1)          # velocidad normal en aristas, paredes = 0
 
    def initial_bump(self, x0=0.3, sigma=0.05, amp=0.05):
        x = (np.arange(self.N) + 0.5) * self.dx
        self.eta = amp * np.exp(-((x - x0) ** 2) / (2 * sigma ** 2))
        self.u[:] = 0.0
 
    def step(self, dt):
        N, dx, g, H, th = self.N, self.dx, self.g, self.H, self.theta
        k = g * H * dt * dt * th * th / (dx * dx)   # coeficiente de la ecuación de onda
        alpha = H * dt / dx
 
        # G_j: parte explícita reunida en las aristas (paredes = 0)
        G = np.zeros(N + 1)
        G[1:N] = self.u[1:N] - g * (dt / dx) * (1 - th) * (self.eta[1:] - self.eta[:-1])
 
        # Tridiagonal: -k η_{i-1} + (1+2k) η_i - k η_{i+1} = R_i
        ab = np.zeros((3, N))
        ab[0, 1:] = -k                         # supra-diagonal
        ab[2, :-1] = -k                        # sub-diagonal
        ab[1, :] = 1.0 + 2 * k
        ab[1, 0] = 1.0 + k                     # pared izquierda
        ab[1, -1] = 1.0 + k                    # pared derecha
 
        rhs = self.eta.copy()
        rhs -= alpha * th * (G[1:] - G[:-1])
        rhs -= alpha * (1 - th) * (self.u[1:] - self.u[:-1])
 
        new_eta = solve_banded((1, 1), ab, rhs)
 
        new_u = self.u.copy()
        new_u[1:N] = self.u[1:N] - g * (dt / dx) * (
            th * (new_eta[1:] - new_eta[:-1]) + (1 - th) * (self.eta[1:] - self.eta[:-1])
        )
        new_u[0] = new_u[-1] = 0.0
        self.eta, self.u = new_eta, new_u
        return self.eta
 
if __name__ == "__main__":
    sim = CasulliShallowWater1D(N=200, theta=0.6)
    sim.initial_bump()
    c = np.sqrt(sim.g * sim.H)
    dt_explicit_limit = sim.dx / c
    dt = 4.0 * dt_explicit_limit                # 4× el límite explícito
    for _ in range(int(0.6 / dt)):
        sim.step(dt)
    print(f"max|η| = {np.max(np.abs(sim.eta)):.4f}  (sigue acotado)")

Un paso explícito diverge en $\Delta t = 4 \Delta t_\text{ex}$ . Casulli termina el intervalo de longitud 0,6 con ese mismo paso. La fase se atrasa. La amplitud muere más rápido cuando θ se acerca a 1. θ=0,5 parte los pesos por la mitad y es de segundo orden. Pasado 0,6 el esquema es de primer orden. El valor por defecto del código, 0,6, ya está de ese lado de primer orden.

Con un paso de tiempo mayor#

Longitud 1 y 200 celdas dan $\Delta x = 0{,}005$ . Con $c=\sqrt{9{,}81}\approx 3{,}13$ m/s el límite explícito es unos 0,0016 s, y cuatro veces eso es unos 0,0064 s. Al descartar la fracción de 0,6 dividido entre ese paso quedan 93 pasos. La salida solo dice que max|η| siguió finito. No dice cuánto decayó el montículo inicial de 0,05.

Ese 0,05 es el 5% de la profundidad $H=1$ . La linealización $\partial_t \eta + H\partial_x u = 0$ se cierra solo mientras $\eta$ queda bien por debajo de $H$ . La celeridad es entonces un solo número, $\sqrt{gH}$ , y la $k$ tridiagonal es la misma en cada celda. Si el montículo llega al orden de la profundidad, $c$ cambia de un sitio a otro. Esta $k$ constante ya no cierra la ecuación de onda.

Con θ=0,5 el mismo paso cuádruple da $k=16\times 0{,}25=4$ . El excedente diagonal sigue en 1. El borde de estabilidad no es la matriz. Es el peso temporal $(θ-1/2)$ . En cero, la viscosidad numérica se apaga y los modos de malla van y vienen sin decaer.

Lo que falta#

No es un solver universal.

El orden temporal cambia con la perilla. θ=0,5 es de segundo orden, pero está en el borde de estabilidad y los frentes de onda gravitatoria se doblan. Pasado 0,6 el orden es uno. Si la amplitud del mareógrafo es el producto, conviene achicar el paso o pasar a un IMEX de orden mayor.

La advección se queda explícita. El upwind simple es pegajoso. Por eso un trazado Euleriano–Lagrangiano recorre la trayectoria hacia atrás e interpola. El paso puede cruzar más de una celda y aun así sostenerse, pero el resultado se mueve con el orden de interpolación y con el error del punto de partida.

El balance hidrostático $\partial p/\partial z = -\rho g$ está supuesto. En rompientes y caídas, donde la aceleración vertical es grande, hace falta una corrección no hidrostática.

Un triángulo obtuso y el paso de cinco minutos se retira#

Cinco minutos con θ=0,6 compran un recorte de 1 500 veces en el número de pasos y pagan un amortiguamiento de primer orden de la altura de onda. Si la amplitud del registro de marea es el producto, no conviene empujar θ hacia 1. El excedente diagonal ya es 1, así que un paso más implícito no vuelve la matriz más segura. Lo que muere es la onda.

La diferencia iguala $\partial \eta / \partial n$ solo mientras la línea de centros es perpendicular a la cara. Un triángulo obtuso empuja el circuncentro fuera de la celda y esa identidad termina. Queda una pieza tangencial. Si esa pieza puede reponerse como se repone un flujo difusivo es la pregunta de la corrección del flujo difusivo en una cara no ortogonal.