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

Celdas mojadas y secas en un mismo sistema, con convergencia monótona asegurada.

$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.

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.

\mathbf{V}(\eta) + T\eta = \mathbf{b}

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).

a_i(\eta) = p_i(\eta) - q_i(\eta), \quad p_i, q_i \ge 0,\ \text{no decrecientes}

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

V_i(\eta) = V_{1,i}(\eta) - V_{2,i}(\eta), \quad V_{j,i}(\eta) = \int_{-\infty}^\eta p_i\ \text{o}\ q_i

$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:

V_1(\eta) - V_2(\eta) + T\eta = \mathbf{b}

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

V_1(\eta^n) + (T - Q^{n-1})\eta^n = \mathbf{b} + V_2^{n-1} - Q^{n-1}\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$ :

(T + P^{n,m-1} - Q^{n-1})\eta^{n,m} = P^{n,m-1}\eta^{n,m-1} - V_1^{n,m-1} + \mathbf{d}^{n-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". Vale la pena moverlo con la mano:

T = 0.50b = 0.90

step 1/4 · outer n=0 · eta=-0.4000 · |r|=1.10e+0

Pink = outer iterate (rises from below); green = inner iterate (falls from above). They squeeze toward the true root of V(eta)+T*eta=b with V(eta)=max(0,min(1,eta)).

Si T baja a 0.05, $T+P-Q$ se vuelve casi singular y las iteradas igual terminan en cinco pasos. Si b se empuja más allá de 1.5, la raíz cae en la cola plana de $V$ (donde Newton oscilaría sin fin) — el bucle interior frena en un paso.

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

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 mismo algoritmo, ejecutado día a día sobre un corte 1D del acuífero. La fuerza de la bomba es un control deslizante; el botón "+1 day" avanza un día.

pump flux phi = 0.0008 (m/s per cell)

day = 0 · pressurized cells = 58 · phreatic = 0 · dry = 2 · total water V = 3998.9 m

Cyan = pressurized (confined) region; teal = phreatic; black = dry. Yellow line is piezometric head eta(x). Pump near the center steadily drains the aquifer; the wet/dry boundary advances inward as days pass.

Subiendo la bomba por encima de 0.0025 y dejándolo correr diez días, un parche negro de tierra seca emerge en el centro y se expande. El solver es indiferente al cambio diario del conjunto de celdas mojadas. Con Newton clásico habría que reconstruir el Jacobiano cada vez que el estencil cambia. El bucle anidado resuelve sin retocar una sola línea.

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

No es una panacea. Algunas trampas:

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.
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.
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í.
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.

Próximas lecturas#

Brugnano & Casulli (2008/2009) — teoremas de terminación en pasos finitos para sistemas lineales por tramos. Ancestro directo del presente trabajo.
Casulli (2009) — algoritmo de wetting/drying de alta resolución del mismo autor.
Casulli & Zanolli (2010) — Newton anidado aplicado a Richards en forma mixta (caso $a(\eta)$ suave).