Frontal polymerization Glass Fiber 2d Hires

def _kamal_sourour_rate_vec(T: np.ndarray, alpha: np.ndarray) -> np.ndarray:
    """Vectorised Kamal-Sourour autocatalytic cure rate.

    Clamps inputs for numerical safety before evaluating the Arrhenius
    terms.  Returns dα/dt for each node.
    """
    T_safe = np.clip(T, _T_FLOOR, _T_CEIL)
    alpha_safe = np.clip(alpha, 0.0, 1.0)

    inv_RT = 1.0 / (_R_GAS * T_safe)
    arg1 = np.clip(_E1 * inv_RT, 0.0, _EXP_ARG_MAX)
    arg2 = np.clip(_E2 * inv_RT, 0.0, _EXP_ARG_MAX)

    k1 = _A1 * np.exp(-arg1)
    k2 = _A2 * np.exp(-arg2)

    return (k1 + k2 * np.power(alpha_safe, _M)) * np.power(1.0 - alpha_safe, _N_ORD)

def _fp_glass_fiber_2d_hires_rhs(t, y):
    T_flat = np.clip(y[:_GF_N2D], _T_FLOOR, _T_CEIL)
    alpha_flat = np.clip(y[_GF_N2D:], 0.0, 1.0)

    T = T_flat.reshape(_GF_NX, _GF_NY)
    alpha = alpha_flat.reshape(_GF_NX, _GF_NY)

    dadt_2d = _kamal_sourour_rate_vec(T, alpha)

    dT = np.empty((_GF_NX, _GF_NY))

    for i in range(_GF_NX):
        for j in range(_GF_NY):
            # --- x-direction (fiber) Laplacian ---
            if j == 0:
                # Left edge: Dirichlet T = T_LEFT (held constant)
                lap_x = 0.0
            elif j == _GF_NY - 1:
                # Right edge: convective BC
                # Ghost: T_ghost = T[i,j] + (h*dx/k_x)*(T_amb - T[i,j])
                T_ghost = T[i, j] + (_H_CONV * _GF_DX / _K_FIBER_G) * (_T_AMBIENT - T[i, j])
                lap_x = (T[i, j - 1] - 2.0 * T[i, j] + T_ghost) * _GF_INV_DX2
            else:
                T_left = _GF_T_LEFT if j == 1 else T[i, j - 1]
                lap_x = (T_left - 2.0 * T[i, j] + T[i, j + 1]) * _GF_INV_DX2

            # --- y-direction (transverse) Laplacian ---
            if i == 0:
                # Top edge (row 0): convective BC
                T_ghost = T[i, j] + (_H_CONV * _GF_DY / _K_TRANS_G) * (_T_AMBIENT - T[i, j])
                lap_y = (T_ghost - 2.0 * T[i, j] + T[i + 1, j]) * _GF_INV_DY2
            elif i == _GF_NX - 1:
                # Bottom edge: convective BC
                T_ghost = T[i, j] + (_H_CONV * _GF_DY / _K_TRANS_G) * (_T_AMBIENT - T[i, j])
                lap_y = (T[i - 1, j] - 2.0 * T[i, j] + T_ghost) * _GF_INV_DY2
            else:
                lap_y = (T[i - 1, j] - 2.0 * T[i, j] + T[i + 1, j]) * _GF_INV_DY2

            dT[i, j] = _GF_DIFF_X * lap_x + _GF_DIFF_Y * lap_y + _SRC_COEFF * dadt_2d[i, j]

    # Left edge (j=0) is Dirichlet — temperature held constant
    dT[:, 0] = 0.0

    dy = np.empty(_GF_DIM)
    dy[:_GF_N2D] = dT.ravel()
    dy[_GF_N2D:] = dadt_2d.ravel()
    return dy