# Verificación de parada y reinicio de la mezcla (revisión A1 §2–§3).
# Plano inclinado uniforme de pendiente S. Teoría (flujo de Bingham integrado en profundidad, régimen de fluencia):
# una capa en reposo es estable si ρ g h S_superficie ≤ τ_y  ⇒  espesor crítico h_c = τ_y / (ρ g S).
# Caso 1 (parada): se libera un bloque de 2 m de espesor; al detenerse, el espesor en la parte media del depósito debe tender a h_c
#   (en la cola, más delgado; en el frente, más grueso por la pendiente de superficie). Se reporta además el transporte residual.
# Caso 2 (reinicio): a t = 900 s se agrega mezcla aguas arriba; el depósito debe volver a moverse y detenerse de nuevo.
# Caso 3 (convergencia temporal del operador): CFL 1,0 vs 0,5 vs 0,25 → alcance del frente.
import sys, json, math, numpy as np, anuga
RHO, G, S, TAU, ETA, K, N = 1700.0, 9.81, 0.05, 100.0, 1.0, 100.0, 0.025
HC = TAU / (RHO * G * S)

def corrida(cfl=1.0, reinicio=True, L=600.0, W=20.0, dx=2.0):
    dom = anuga.rectangular_cross_domain(int(L / dx), int(W / dx), len1=L, len2=W)
    dom.set_flow_algorithm('DE1'); dom.set_CFL(cfl); dom.set_minimum_storable_height(1e-3)
    dom.set_quantity('elevation', lambda x, y: -S * x); dom.set_quantity('friction', 0.0)
    dom.set_quantity('stage', lambda x, y: -S * x + np.where((x > 20) & (x < 40), 2.0, 0.0))
    Br = anuga.Reflective_boundary(dom); dom.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
    Qf = (lambda t: 2.0 if 900 <= t <= 960 else 0.0) if reinicio else (lambda t: 0.0)
    anuga.Inlet_operator(dom, anuga.Region(dom, center=[10.0, W / 2], radius=4.0), Q=Qf, default=0.0)
    class Reo(anuga.Operator):
        def __call__(self):
            dt = self.domain.get_timestep(); h = np.maximum(self.stage_c - self.elev_c, 0); wet = h > 1e-3; hs = np.where(wet, h, 1.0)
            hu, hv = self.xmom_c, self.ymom_c; Us = np.where(wet, np.hypot(hu, hv) / hs, 0)
            A_ = TAU / (RHO * hs); B_ = K * ETA / (8 * RHO * hs ** 2); C_ = G * N ** 2 / hs ** (4 / 3)
            r = np.maximum(0, Us - dt * A_); Un = 2 * r / (1 + B_ * dt + np.sqrt((1 + B_ * dt) ** 2 + 4 * C_ * dt * r))
            f = np.where(wet & (Us > 0), Un / np.where(Us > 0, Us, 1), 0); hu *= f; hv *= f
        def parallel_safe(self): return True
        def statistics(self): return ''
        def timestepping_statistics(self): return ''
    Reo(dom)
    x = dom.centroid_coordinates[:, 0]; hist = []
    for t in dom.evolve(yieldstep=30.0, finaltime=1800.0):
        h = np.maximum(dom.quantities['stage'].centroid_values - dom.quantities['elevation'].centroid_values, 0)
        mom = np.hypot(dom.quantities['xmomentum'].centroid_values, dom.quantities['ymomentum'].centroid_values)
        front = float(x[h > 0.01].max()) if (h > 0.01).any() else 0.0
        hist.append((float(t), front, float(np.sum(h * dom.areas)), float(np.sum(mom * dom.areas))))
    h = np.maximum(dom.quantities['stage'].centroid_values - dom.quantities['elevation'].centroid_values, 0)
    dep = x[h > 0.01]; mid = (h > 0.01) & (x > np.percentile(dep, 35)) & (x < np.percentile(dep, 65))
    return np.array(hist), float(np.median(h[mid])) if mid.any() else 0.0

if __name__ == '__main__':
    out = {'h_c_teorico_m': round(HC, 4), 'parametros': dict(rho=RHO, S=S, tau_y=TAU, eta=ETA, K=K, n=N)}
    hist, hmid = corrida(reinicio=True)
    i880 = np.searchsorted(hist[:, 0], 880); i1800 = len(hist) - 1
    out['parada'] = {'frente_m_t0_300_880': [round(hist[np.searchsorted(hist[:, 0], tt), 1], 2) for tt in (0, 300, 880)],
                     'deriva_frente_600_a_880_m': round(hist[i880, 1] - hist[np.searchsorted(hist[:, 0], 600), 1], 3),
                     'momento_total_880': round(hist[i880, 3], 4)}
    out['reinicio'] = {'frente_antes_m': round(hist[i880, 1], 2), 'frente_final_m': round(hist[i1800, 1], 2),
                       'se_movio': bool(hist[i1800, 1] > hist[i880, 1] + 1.0),
                       'volumen_final_m3': round(hist[i1800, 2], 2), 'volumen_esperado_m3': round(2.0 * 20 * 20 + 2.0 * 60, 2)}
    out['espesor_medio_deposito_final_m'] = round(hmid, 4)
    out['convergencia_CFL'] = {}
    for cfl in (1.0, 0.5, 0.25):
        hh, hm = corrida(cfl=cfl, reinicio=False)
        out['convergencia_CFL'][str(cfl)] = {'frente_final_m': round(hh[-1, 1], 2), 'espesor_medio_m': round(hm, 4)}
    json.dump(out, open('res/verif_parada.json', 'w'), indent=1, ensure_ascii=False); print(json.dumps(out, indent=1, ensure_ascii=False))