estimation of gradient $$(\mathcal{P}_{\varepsilon}) : \left\{\begin{array}{ll}
\displaystyle -div\left(A(x)\nabla u_\varepsilon(x)\right)= \dfrac{a(x)}{|u_\varepsilon(x)|+\varepsilon} &\mbox{ in }\Omega \\\\
u_\varepsilon = 0 &\mbox{ over }\partial \Omega
\end{array}
\right.$$
Assume that $\Omega \in C^{1,1}$, be a bounded open set,
$a \in L^\infty, a\geq0 \mbox{ and } a\neq 0$ for $x \in \Omega, A(x)=(a_{ij}(x))$
a uniformly coercive matrix with coefficients $a_{i,j} \in C^{0,1}(\bar{\Omega}$). Then, $\mathcal{P}_\varepsilon$ have a unique solution  $u_\varepsilon \in W^1_0\cap W^{1,p}(\Omega)$.
Moroever, $\exists\  C_\Omega$ uniformly independent of $\varepsilon$ such that
\begin{equation}
\| \nabla u_\varepsilon \|_{L^2} \leqslant C_\Omega \int_\Omega a(x)\,dx.
\end{equation}
I'm not sure how to get this inequality, is there any help?
 A: I am a bit rusty on this, so work with me. We are not proving existence or uniqueness, only answering from where the inequality might come. Since $A(x)$ is uniformly coercive, that means that since
$${\rm -div}(A(x)\nabla u_\epsilon)\ge0,$$ then we can apply  maximum principle to $u_\epsilon$ (is super-harmonic). Therefore, the minimum of $u_\epsilon$ is attained in $\partial\Omega$, where it is zero.
Therefore $u_\epsilon>0$ in $\Omega$ and $|u_\epsilon|=u_\epsilon$. We can re-write the original problem as
$$ -(u_\epsilon+\epsilon)\,{\rm div}(A\nabla u_\epsilon)=a~~.\quad(1)$$
Since 
$$f\, {\rm div}(A\nabla f)+(\nabla f)^T A (\nabla f)={\rm div}(fA\nabla f)$$ we have the following generalized Green's identity 
$$\int_\Omega f {\rm div}(A\nabla f)+(\nabla f)^T A (\nabla f)dx=\int_{\partial\Omega}f(A\nabla f)\cdot d\vec{S}.$$
We apply this identity to Eq.(1) and we get, using $u|_{\partial\Omega}=0$,
\begin{align}
 \int_\omega (\nabla u_\epsilon)^T A (\nabla u_\epsilon)-\epsilon{\rm div}(A(x)\nabla u_\epsilon)dx&=\int_\Omega a(x)dx\\
 \int_\omega (\nabla u_\epsilon)^T A (\nabla u_\epsilon)&=\int_\Omega a(x)\left(1-\frac{\epsilon}{\epsilon+u_\epsilon}\right)dx\\&\le\int_\Omega a(x)dx~~, 
\end{align}
where in the second step we use that $u_\epsilon$ is solution. But $\alpha|v|^2\le|v^TAv|\le \beta |v|^2$ for two constants $0<\alpha\le\beta$ in general since $A$ components are continuous in $\bar{\Omega}$ and form coercive matrix. Therefore $$ ||\nabla u_\epsilon||^2_2\le\frac{1}{\alpha}\int_\Omega a(x)dx, $$ where constant $\alpha$ is uniform.
