  For time-dependent PDEs, the numerical schemes can be rendered
bound-preserving without losing conservation and accuracy, by a post processing
procedure of solving a constrained minimization in each time step. Such a
constrained optimization can be formulated as a nonsmooth convex minimization,
which can be efficiently solved by first order optimization methods, if using
the optimal algorithm parameters. By analyzing the asymptotic linear
convergence rate of the generalized Douglas-Rachford splitting method, optimal
algorithm parameters can be approximately expressed as a simple function of the
number of out-of-bounds cells. We demonstrate the efficiency of this simple
choice of algorithm parameters by applying such a limiter to cell averages of a
discontinuous Galerkin scheme solving phase field equations for 3D demanding
problems. Numerical tests on a sophisticated 3D Cahn-Hilliard-Navier-Stokes
system indicate that the limiter is high order accurate, very efficient, and
well-suited for large-scale simulations. For each time step, it takes at most
$20$ iterations for the Douglas-Rachford splitting to enforce bounds and
conservation up to the round-off error, for which the computational cost is at
most $80N$ with $N$ being the total number of cells.
