This is a nice question. In my view, the choice of test function is motivated by the idea of taking some positive quantity that is a supersolution to an elliptic equation, and looking at the equation for its logarithm. The new equation contains a useful term that is quadratic in the gradient. This idea is pervasive in geometry and elliptic PDEs, and some examples are below.
(1) The basic case to consider is that $u$ is positive and superharmonic. Then $v := -\log u$ satisfies $|\nabla v|^2 \leq \Delta v$, which gives local universal bounds on the integral of $|\nabla v|^2$ after multiplying by standard cutoffs and integrating by parts.
(2) In your context, the choice of test function $H(u)$ satisfies $$a^{ij}\partial_iu\partial_j(H(u)) = a^{ij}\partial_i(V(u))\partial_j(V(u)),$$ where $V(u) = c_1\log(c_2 - u)$ with $c_2 - u$ positive. I view the estimate as coming from integrating the equation for $V(u)$. To illustrate how this works in a simple context, assume that $u \in C^2_0(B_1)$ satisfies $\Delta u \geq -1$, and that $u \leq M$. Then for $w := M+1-u$ we have that $v := -\log(w)$ satisfies $|\nabla v|^2 \leq 1 + \Delta v$. Thus the integral of $|\nabla v|^2$ is bounded by $1$.
(3) The Bombieri-De Giorgi-Miranda interior gradient estimate for a solution $u$ to the minimal surface equation is based on the fact that the vertical component $\nu^{n+1}$ of the unit normal to the graph of $u$ is positive and superharmonic (on the graph). The proof uses the equation for $v := -\log(\nu^{n+1})$, which just as above contains a useful term quadratic in $|\nabla v|$.
(4) The Li-Yau proof of the Harnack inequality for a harmonic function $u$ is obtained by looking at the quantity $w := |\nabla (-\log u)|^2$. The key is that $w$ solves a differential inequality with the powerful term $\frac{2}{n}w^2$, which allows one to bound $w$ from above locally by a universal constant.
I am sure there are many other interesting examples, and I am not sure where the first instances of the "log trick" appeared.