FP Kamal-Sourour 1D Front (dim=20)

def _kamal_sourour_rate(T: float, alpha: float) -> float:
    """Kamal-Sourour autocatalytic cure rate at a single point.

    Guards against unphysical states: alpha is clamped to [0, 1] and T is
    floored at _T_FLOOR to prevent Arrhenius overflow from negative or
    near-zero temperatures during stiff integration.
    """
    alpha = min(max(alpha, 0.0), 1.0)
    T = max(T, _T_FLOOR)
    inv_RT = 1.0 / (_R_GAS * T)
    k1 = _A1 * np.exp(-_E1 * inv_RT)
    k2 = _A2 * np.exp(-_E2 * inv_RT)
    return (k1 + k2 * alpha ** _M) * (1.0 - alpha) ** _N

def fp_kamal_sourour_1d_rhs(t, y):
    T = y[:_N1D].copy()
    alpha = y[_N1D:].copy()

    np.clip(T, _T_FLOOR, None, out=T)
    np.clip(alpha, 0.0, 1.0, out=alpha)

    dT = np.empty(_N1D)
    dalpha = np.empty(_N1D)

    for i in range(_N1D):
        R_i = _kamal_sourour_rate(T[i], alpha[i])
        dalpha[i] = R_i

        # Laplacian with Dirichlet BCs at both ends
        T_left = _T_IGNITION if i == 0 else T[i - 1]
        T_right = _T_AMBIENT if i == _N1D - 1 else T[i + 1]
        lap = (T_left - 2.0 * T[i] + T_right) * _INV_DX2_1D

        dT[i] = _DIFF_1D * lap + _SRC_COEFF_1D * R_i

    dy = np.empty(2 * _N1D)
    dy[:_N1D] = dT
    dy[_N1D:] = dalpha
    return dy