Poroviscoelasticity involves primary and secondary consolidation of fluid-filled porous media, i.e., the superposition of instantaneous and delayed responses to loadings. In this talk, we cast the corresponding model as a gradient flow depicting the dissipative character of the system. By exploiting the inherent minimization structure at semidiscrete level, we develop and analyze novel robust block-partitioned solvers extending the popular fixed-stress and undrained splits (for poroelasticity). The systematic approach demonstrated here is applicable to a large class of coupled non-linear models.