FP Thermal Front 2D Composite (dim=50)

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_thermal_front_2d_rhs(t, y):
    T_flat = y[:_N2D].copy()
    alpha_flat = y[_N2D:].copy()

    np.clip(T_flat, _T_FLOOR, None, out=T_flat)
    np.clip(alpha_flat, 0.0, 1.0, out=alpha_flat)

    T = T_flat.reshape(_NX, _NY)
    alpha = alpha_flat.reshape(_NX, _NY)

    dT = np.empty((_NX, _NY))
    dalpha = np.empty((_NX, _NY))

    for i in range(_NX):
        for j in range(_NY):
            R_ij = _kamal_sourour_rate(T[i, j], alpha[i, j])
            dalpha[i, j] = R_ij

            # --- x-direction (fiber) Laplacian ---
            if i == 0:
                # Left edge: Dirichlet T = T_IGNITION (held constant, so
                # dT/dt = 0 but the Laplacian stencil still uses this value
                # as a ghost).  The actual derivative for i=0 nodes is forced
                # to zero below.
                lap_x = 0.0
            elif i == _NX - 1:
                # Right edge: Neumann dT/dx=0 -> ghost T[NX, j] = T[NX-1, j]
                lap_x = (T[i - 1, j] - T[i, j]) * _INV_DX2_2D
            else:
                lap_x = (T[i - 1, j] - 2.0 * T[i, j] + T[i + 1, j]) * _INV_DX2_2D

            # --- y-direction (transverse) Laplacian ---
            if j == 0:
                # Bottom edge: Neumann dT/dy=0 -> ghost T[i, -1] = T[i, 0]
                lap_y = (T[i, j + 1] - T[i, j]) * _INV_DY2_2D
            elif j == _NY - 1:
                # Top edge: Neumann dT/dy=0 -> ghost T[i, NY] = T[i, NY-1]
                lap_y = (T[i, j - 1] - T[i, j]) * _INV_DY2_2D
            else:
                lap_y = (T[i, j - 1] - 2.0 * T[i, j] + T[i, j + 1]) * _INV_DY2_2D

            dT[i, j] = _DIFF_X * lap_x + _DIFF_Y * lap_y + _SRC_COEFF_2D * R_ij

    # Left edge is Dirichlet — temperature is held constant by the BC
    dT[0, :] = 0.0

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