We work on a closed Riemannian manifold $M$. Let $u$ and $v$ be the non-negative weak solutions of $$au_t - 2d\,\Delta au = cv - f$$ $$bv_t - d\,\Delta bv = f$$ $$u(0)=u_0, \quad v(0)=v_0$$ where $f$ is a given non-negative function in $L^2(0,T;H^1(M))$ and $a$, $b$ , $c$ and $d$ are positive constants, and $u_0$ and $v_0$ are given non-negative functions that are bounded on $M$.

These solutions lie in the space $L^2(0,T;H^1(M))$ with time derivative in $L^2(0,T;L^2(M))$.

Are $u$ and $v$ in $L^\infty(0,T;L^\infty(M))$?

I'd like to prove this using weak/variational techniques and not recourse to more "classical" theory if possible.

Note that $f$ is not bounded. The problem is that the diffusion coefficients are not equal. If they were, we can add the two equations and test with $(au+bv-C)^+$ where $C$ depends on the $L^\infty$ bounds on the initial data. (The term $cv$ on the right hand side of the first equation can be neglected when we add the equations because we can rescale $v$ by an expontential function and so get rid of it).