[Revisión de artículo] Partir Newton en dos y la superficie libre se vuelve resoluble — Casulli–Zanolli (2012)

$V(\eta) + T\eta = \mathbf{b}$. Una sola línea. La primera vez parece que basta un paso de Newton para cerrarla. Hasta que aparece el borde del pozo — celdas cuya profundidad cruza el cero. Allí el Jacobiano se anula sin avisar. La misma celda mojada en un instante y seca en el siguiente (wetting/drying). Es donde tropieza cualquier código de superficie libre basado en Newton. Lo que Casulli y Zanolli hicieron en 2012 fue partir Newton en dos. La iteración exterior sube desde abajo, la interior baja desde arriba, y aprietan la solución entre las dos. La convergencia monótona queda garantizada.

Conviene bajar T hasta 0.05 y empujar b más allá de 1.5. El Newton anidado escalar aprieta desde arriba y desde abajo a la vez.

Con T en 0.05 la matriz $T+P-Q$ queda casi singular y la iteración igual se detiene en cinco pasos. Pasado b = 1.5 la raíz se sienta en la cola plana de $V$, donde $\eta > 1$, y el lazo interior frena en un paso.

Datos del artículo#

• Autores: Vincenzo Casulli, Paola Zanolli
• Revista: Journal of Computational and Applied Mathematics, vol. 236 (2012), pp. 3937–3947
• DOI: 10.1016/j.cam.2012.04.025
• Palabras clave: mildly nonlinear systems / wetting and drying / free-surface hydrodynamics / Richards equation / confined–unconfined aquifers

Q1. ¿Por qué Newton muere tanto en superficies libres?#

La ecuación es de aspecto inofensivo.

Aquí $\eta \in \mathbb{R}^N$ es la altura piezométrica por celda, $T$ es una matriz simétrica semidefinida positiva (estencil de difusión discreto) y $\mathbf{V}$ es una función no lineal diagonal — el volumen por celda $V_i(\eta_i) = \int_{-\infty}^{\eta_i} a_i(z)\,dz$. La función $a_i$ solo necesita ser de variación acotada (se permiten discontinuidades). Tres ejemplos canónicos:

• wetting/drying: $V(\eta) = \max(0, \eta)$. Las celdas secas tienen $V'=0$.
• acuífero confinado–no confinado: $V(\eta) = \max(0, \min(c-(-h), \eta-(-h)))$. Dos quiebres.
• ecuación de Richards en forma mixta: $V(\eta)$ depende del suelo, suave pero casi plana en intervalos amplios.

Los tres comparten un rasgo: $V'(\eta)$ se anula sobre intervalos enteros. El Newton estándar utiliza $J = V'(\eta^{(k)}) + T$. Cuando $V'=0$ y la fila correspondiente de $T$ acopla débilmente (o roza el núcleo de $T$), $J$ se vuelve singular. Si el valor inicial no está suficientemente cerca de la verdad, la iteración diverge — esa es la plaga de tres décadas de los códigos de superficie libre.

Q2. Cortar $\mathbf{V}(\eta)$ en dos — descomposición de Jordan#

El primer truco es descomponer la función. Variación acotada de $a_i$ significa que se escribe como diferencia de dos funciones no decrecientes (descomposición de Jordan).

Automáticamente $V_i$ también se descompone.

$V_1$ y $V_2$ son ambas no decrecientes y acotadas. En el caso wetting/drying: $p(\eta)=1\,[\eta \ge 0]$, $q\equiv 0$. En el caso confinado–no confinado: $p(\eta)=1\,[\eta \ge -h]$, $q(\eta)=1\,[\eta > c]$. Una vez separada, cada pieza es monótona — esa frase es la semilla de todo lo demás.

Q3. Dos linealizaciones, una dentro y otra fuera#

El sistema a resolver:

El bucle exterior lineariza $V_2$ alrededor de $\eta^{n-1}$.

con $Q^{n-1} = \mathrm{diag}(q_i(\eta_i^{n-1}))$. Punto de partida: $\eta^0 \le \ell$ (valla inferior). Esto sigue siendo no lineal, así que el bucle interior lineariza $V_1$:

con $\eta^{n,0} \ge u$ (valla superior). Las vallas $\ell, u$ se escogen para que $T+P-Q$ resulte una matriz de Stieltjes (M-matriz simétrica). El premio:

• Las iteradas internas son monótonas decrecientes desde arriba ($\eta^{n,m+1} \le \eta^{n,m}$).
• Las iteradas externas son monótonas crecientes desde abajo ($\eta^{n+1} \ge \eta^n$).
• Se aprietan sobre la única solución.

La convergencia se demuestra con desigualdades, no con argumentos de aritmética flotante. Desaparece el viejo defecto "Newton converge sólo cerca de la verdad". Los cinco pasos y el paso único del control de arriba son esas desigualdades, contadas. Si la matriz interior sigue siendo de Stieltjes, el sistema es el mismo problema SPD que un ciclo V que limpia las ondas largas.

Q4. Un pozo con bomba en 1D — implementación en Python#

La §6 del artículo estudia un acuífero confinado–no confinado en 2D. Lo reducimos a 1D y lo escribimos desde cero. El problema de prueba es fondo/techo paraboloidal más una bomba al centro: con el correr de los días la región freática crece y luego aparece la zona seca. Ninguno de los problemas estándar (Sod, Karman, blast) encaja, así que se conserva la prueba del propio artículo.

import numpy as np
 
# Volumen por celda: V_i(eta) = a0 * clamp(eta, -h_i, c_i) (trasladado por -h_i)
def cell_volume(eta, h, c, a0=0.3):
    return a0 * (np.clip(eta, -h, c) - (-h))
 
# Jacobianos diagonales del corte de Jordan
def p_diag(eta, h, a0=0.3):
    # derivada de V1 = a0 * max(0, eta - (-h))
    return np.where(eta >= -h, a0, 0.0)
 
def q_diag(eta, c, a0=0.3):
    # derivada de V2 = a0 * max(0, eta - c)
    return np.where(eta > c, a0, 0.0)
 
def nested_newton_step(eta_prev, h, c, T, dt, phi, tol=1e-6):
    """Un paso temporal con Newton anidado monótono.
 
    Resuelve V(eta_new) + T eta_new = V(eta_prev) + dt*phi,
    donde T es la matriz SPD de -dt*Laplaciano.
    """
    N = len(eta_prev)
    V_prev = cell_volume(eta_prev, h, c)
    rhs0 = V_prev + dt * phi
 
    # Arranque exterior: muy por debajo de la valla baja (todo seco)
    eta = np.full(N, -h.max() - 1.0)
    for n_out in range(20):
        Q = np.diag(q_diag(eta, c))
        V2 = np.where(eta > c, 0.3 * (eta - c), 0.0)
        d = rhs0 + V2 - Q @ eta
 
        # Arranque interior: muy por encima de la valla alta (todo presurizado)
        eta_in = np.full(N, c.max() + 1.0)
        for m_in in range(30):
            P = np.diag(p_diag(eta_in, h))
            V1 = np.where(eta_in > -h, 0.3 * (eta_in - (-h)), 0.0)
            f = P @ eta_in - V1 + d
            A = T + P - Q                     # Stieltjes (M-matriz simetrica)
            eta_in = np.linalg.solve(A, f)    # 1D denso; 2D requiere PCG
            r_in = V1 + (T - Q) @ eta_in - d
            if np.max(np.abs(r_in)) < tol:
                break
        eta = eta_in
        V_cur = cell_volume(eta, h, c)
        r_out = V_cur + T @ eta - rhs0
        if np.max(np.abs(r_out)) < tol:
            return eta, n_out + 1, m_in + 1
    return eta, 20, m_in + 1
 
# Pozo 1D (corte paraboloidal)
N = 41
L = 1000.0
x = np.linspace(0, L, N)
h_arr = 10 * np.maximum(0, 1 - ((x - L/2) / (L/2))**2)
c_arr = h_arr.copy()
 
# Matriz de difusion T = dt*kappa/dx^2 * Laplaciano (flujo nulo en bordes)
dx = L / (N - 1)
dt = 86400.0     # 1 dia
kappa = 1.0
alpha = dt * kappa / dx**2
T = -alpha * np.eye(N, k=-1) + 2*alpha*np.eye(N) - alpha*np.eye(N, k=1)
T[0, 0] = alpha; T[0, 1] = -alpha
T[-1, -1] = alpha; T[-1, -2] = -alpha
 
# Estado inicial: presurizado hasta el techo
eta = c_arr.copy()
phi = np.zeros(N); phi[N//2] = -8e-4   # bomba central (sumidero)
 
for day in range(10):
    eta, n_out, m_in = nested_newton_step(eta, h_arr, c_arr, T, dt, phi)
    pre = int(np.sum(eta >= c_arr - 1e-6))
    dry = int(np.sum(eta <= -h_arr + 1e-6))
    print(f"day {day+1}: pre={pre} dry={dry} outer={n_out} inner_last={m_in}")

Una corrida típica imprime:

day 1: pre=39 dry=0 outer=5 inner_last=5
day 2: pre=37 dry=0 outer=4 inner_last=4
...
day 9: pre=0 dry=8 outer=1 inner_last=5
day 10: pre=0 dry=12 outer=1 inner_last=5

Al subir la bomba y avanzar con +1 day, los conteos pre y dry del listado se ven como una mancha en el corte.

Por encima de 0.0025, diez días abren un parche negro seco en el centro. El número de celdas activas cambia cada día. El lazo anidado no rearma el Jacobiano por eso.

El bucle exterior promedia 3 iteraciones, el interior unas 5. Coincide cualitativamente con la Tabla 1 del artículo. Lo más llamativo es la transición presurizado → freático: ambas iteraciones detienen a menudo en un solo paso (Remark 6 del paper: si $u \le \tilde\eta^{1,1} \le \ell$, la raíz ya está encontrada).

Q5. Evolución del corte 2D paso a paso#

El registro de diez días se lee solo por conteos. Día 1: 39 celdas presurizadas de 41, 0 secas, exterior 5, último interior 5. Día 2: las presurizadas bajan a 37 y el exterior a 4. Los días que el listado omite no se estiman. Día 9: presurizadas 0, secas 8, exterior 1, interior 5. Día 10: solo las secas suben a 12 y el exterior sigue en 1. La bomba es $\phi=-8\times 10^{-4}$ en la celda central, el paso es 86400 y la longitud 1000 se parte en 41 celdas. $a_0=0.3$ entra en $P$ solo en el tramo mojado.

Mientras la carga toca el techo $c$, $q$ está activo y $Q$ vale $a_0$. Bajo $c$, $Q=0$. Bajo el fondo $-h$, $P$ también es 0. La caída del exterior de 5 a 1 en el día 9 acerca el caso al Remark 6. El interior sigue en 5. Las celdas secas añaden filas con $P-Q=0$, y el número de sistemas interiores no baja. Tres pasadas exteriores por unas cinco interiores siguen siendo 15 a 25 sistemas de Stieltjes.

Lo que el algoritmo no puede hacer — un párrafo crítico#

No es una panacea. Algunas trampas:

1. Celdas con $P-Q = 0$. Si una celda seca anula también la fila correspondiente de $T$, la matriz pierde la hipótesis T2 (irreducibilidad). El artículo recomienda excluir las celdas secas del sistema, pero llevar la contabilidad del conjunto activo no es trivial.
2. No es rápido. Por paso temporal hay exterior × interior = 15 a 25 sistemas lineales en promedio. Implicit Euler con CG, cuando converge, suele ser más veloz. El precio que se paga: cero eventos de divergencia.
3. No se extiende directamente al Euler compresible. La técnica vive en el mundo de sistemas levemente no lineales con $V$ diagonal de variación acotada. Las no linealidades de un solver de Riemann no se factorizan así.
4. El solver interior sigue importando. Cada sistema interno es SPD, pero su precondicionamiento es otra disciplina. El paper sugiere PCG; en mallas grandes multigrid suele ser superior.

Tarjeta de reproducibilidad#

• Matemática (§2–§5): autocontenida. Con lápiz y papel se rehacen todas las desigualdades. A.
• Cajas de algoritmo (1 y 2): pseudocódigo completo. A.
• Ejemplo numérico (§6): malla, paso temporal y coeficientes especificados. Treinta líneas de Python bastan para reproducirlo cualitativamente. A.
• Aplicabilidad práctica: OpenFOAM carece de un módulo wetting/drying de primer nivel; los derivados meshFlow/overFlowSolver pueden alojarlo. Los modelos oceánicos y estuáricos (SCHISM, TELEMAC) ya emplean ideas afines.

Cuándo sacar las celdas secas#

Las filas con $P-Q=0$ rompen la hipótesis T2. El artículo pide retirar las celdas secas. El día 10 hay 12 celdas secas de 41. Retirarlas cambia el conjunto activo cada día. Dejarlas mantiene filas nulas. El Newton anidado conserva la convergencia por desigualdades en ambos casos. Lo que cambia es el tamaño del sistema lineal.

La terminación finita de Brugnano y Casulli (2008/2009) vale para sistemas lineales por tramos. Con $a(\eta)$ suave, como en Richards, el lugar es Casulli y Zanolli (2010). Una superficie libre hecha por ondas está en la profundidad que abre una cuña de Kelvin. Aquí la superficie libre es una celda seca.