For this problem: $$ f(x,y) = e^{-x}\, ,0<y<x<\infty $$ Let $Z=2X-Y$, find the density function $f_Z(z)$.
$\underline{Attempt}$
Let $X=\frac{Z+Y}{2}$, then $$ 0<y<x<\infty\implies 0<y<\frac{z+y}{2}<\infty $$ Multiply by $2$ and subtract $y$, obtain $0<y<z$. The mass function is given by $$\begin{aligned} f_Z(z)&=\int_{-\infty}^{\infty}f\left(\frac{z+y}{2},y\right)\; dy \\ &=\int_0^z e^{-\frac{z}{2}-\frac{y}{2}}\; dy\\ &= \left. -2e^{-\frac{z}{2}-\frac{y}{2}}\right\rvert_0^z \\&=2e^{-\frac{z}{2}}-2e^{-z}\; \; (z>0) \end{aligned}$$
However the correct answer is $e^{-\frac{z}{2}}-e^{-z}$. Actually I do not understand what I am doing, or how double integrals work for convolutions at all.
I have a strong foundation in multivariable calculus in calculating area and volumes, but I do not understand how double integrals work for convolutions in probability. For instance, by substituting the variable $x$ for the expression $\frac{z+y}{2}$ and changing the limits of integration, how does it affect the original geometric surface under which volume is supposed to represent the probability?
You were actually very close to the correct solution, while at the same time perhaps quite far because of some misconceptions.
The answer is that the geometry is reflected in the Jacobian, which was missing in your calculation.
Indeed you can consider $Z = 2X - Y = 2X + (-Y)$ as the sum of two dependent random variables $2X$ and $-Y$. Technically, sure, you see the form similar to $\int g(z-y) dy$ so one might call it a convolution in a handwaving manner. However, thinking it as a convolution is misguided since in its basic form convolution works only as the sum of independent random variables.
You were actually making a a 2-dim variable transformation $$\begin{cases} Z = 2X - Y \\ W = Y& \end{cases} \Longleftrightarrow \begin{cases} X = \frac{Z+W}2 & \\ Y = W& \end{cases}$$ where the Jacobian of the inverse transform is $$J = \left|\begin{matrix} \dfrac{ \partial x }{\partial z} & \dfrac{ \partial x }{\partial w} \\ \dfrac{ \partial y }{\partial z} & \dfrac{ \partial y }{\partial w} \end{matrix}\right| = \left|\begin{matrix} \dfrac12 & \dfrac12 \\ 0 & 1 \end{matrix}\right| = \frac{-1}2$$ The 2-dim joint density thus transform as \begin{align} f_{ZW}(z, \, w) &= f_{XY}\bigl(x(z,w), \, y(z,w) \bigr) \cdot |J| \\ &= \frac12 e^{-x(z,w)} \cdot \mathcal{I}_{0<y<x<\infty} \\ &= \frac12 e^{-\frac{z+w}2} \cdot \mathcal{I}_{0 < w < z <\infty} \end{align} where $\mathcal{I}_{\text{region}}$ is the indicator function, basically a shorthand for denoting the domain (support).
Now, we have a 2-dim joint density $f_{ZW}$ so that we can integrate out the auxiliary $W$ to get the marginal $f_Z$.
\begin{align} f_Z(z) &= \int_{w = -\infty}^{\infty} f_{ZW}(z, \, w) \,\mathrm{d}w \\ &= \int_{w = -\infty}^{\infty}\frac12 e^{-\frac{z+w}2} \cdot \mathcal{I}_{0<w<z<\infty} \,\mathrm{d}w \\ &= \frac12 \int_{w = 0}^z e^{-\frac{z+w}2} \,\mathrm{d}w \\ &= \frac12 \left( 2e^{-z/2 } -2e^{-z} \right) \mathcal{I}_{0<z<\infty} \\ &= \left( e^{-z/2 } - e^{-z} \right) \mathcal{I}_{0<z<\infty} \end{align} Technically you were off by just the factor of the absolute value of the Jacobian, $\displaystyle |J| = \left| \frac{-1}2 \right| = \frac12 $. However, using muddled notation (e.g. not setting up a separate notation for the identity mapping $W = Y$) in a 2-dim transform is bad practice at best and plainly incorrect most of the time.
The above variable transform is a standard procedure. Another standard method is directly consider the CDF, as in $\Pr\{2X - Y < z\}$ for a dummy variable $z$, then take the derivative with respect to $z$.
The probability $\Pr\{2X - Y < z\}$ is a 2-dim region integral (over a region of "slanted" triangular region: above $y = 2x -z$ and below $y = x$).