Vertex-patch smoothers are crucial for the robust convergence of geometric multigrid methods with high-order finite elements, yet their adoption is often hindered by the prohibitive computational cost of solving local patch problems, especially in matrix-free settings. In this work, we present a novel vertex patch smoother where local problems are solved inexactly by a nested, matrix-free p-multigrid method. A single iteration of this local solver can be evaluated with O(p^{d+1}) operations, dismantling the traditional trade-off between smoothing effectiveness and computational efficiency. Additionally, we discuss a high-performance implementation that leverages sum-factorization and explicit SIMD vectorization to minimize memory footprint and maximize arithmetic throughput. This implementation achieves optimal O(p^d) memory scaling and rivals the execution speed of simple pointwise smoothers. We further explore the application of this approach to the Stokes equations, demonstrating that a Braess-Sarazin smoother within the p-multigrid framework provides resilience against geometric distortion and large viscosity jumps. All presented methods are based on the concept of fusing residual evaluation with local solving. An open-source finite element library, TensorFEM, inspired by the methods above will be presented. It is a hardware-aware implementation of finite elements: https://github.com/mwichro/TensorFEM