According to the change of variables formula for multivariable calculus, $$d\vec{v}=\left|\det(D\varphi)(\vec{u})\right|d\vec{u}$$ where $\vec{v}=\varphi\vec{u}$ and $\det(D\varphi)(\vec{u})$ is the Jacobian matrix of the partial derivatives of $\varphi$ at the point $\vec{u}$.
How to come up with this relationship (preferably conceptually)?
I'm by no means an expert on this, and this is a question which I've wondered about as well. See also this post, which I think has an excellent answer. I'll rephrase the answer there in my own words, with how I think about it. (btw the treatment of integration I'm following is like that in Spivak's Calculus on manifolds, so when I speak of volume of a set/ measure zero, refer to the definitions provided there)
Suppose $\varphi: A \subset \mathbb{R^n} \to \mathbb{R^n}$ satisfies the hypotheses of the change of variables theorem. Then, the formula \begin{equation} \int_{\varphi(A)} 1 \, dv = \int_A |\text{det}\, D \varphi(u)| \, du \end{equation} expresses how the volume of $\varphi(A)$ is related to "the volume of $A$" (up to a factor of determinant of the derivative). So if we didn't know anything about the change of variables theorem, perhaps the first question we would ask ourselves is "how does the volume of a subset change after we map it by a 'nice' enough $\varphi$". To answer this, note that we have two difficulties, namely the subset $A$ may be complicated, and that the function $\varphi$ may be complicated. A recurring theme in analysis/ calculus is to solve non-linear problems by a well-approximated linear one, and then use a limit argument to show that this is valid. So, we start by simplifying our problem: Suppose $U$ is an open rectangle in $\mathbb{R^n}$, and $\varphi$ is linear. The key insight comes from the fact that in this case \begin{equation} \text{vol }\varphi(U) = |\text{det }\varphi| \cdot \text{vol }U \end{equation}
The proof of this equation is as follows: if $\text{det }\varphi = 0$, then the RHS of the equation is $0$, and the LHS is also $0$, because $\varphi(U)$ has ($n$-dimensional) measure $0$. Otherwise, $\varphi$ can be expressed as a composition of elementary linear transformations, and it is relatively easy to verify that the equation holds for elementary linear transformations (unwinding the definitions of elementary transformation and using Fubini's theorem makes this almost trivial). If it is true for elementary transformations, then it is true for their compositions as well (also easy to verify).
So now in the general case where $A$ and $\varphi$ are arbitrary, the idea is to "chop up" the open set $A$ into a bunch of small rectangles $U$. Then on each of these small rectangles $U$, if $u \in U$ then $\varphi$ is well approximated by its derivative $D\varphi(u)$ so rather than considering the volume of $\varphi(U)$, we can approximate this by the volume of the set $D\varphi(u) [U]$. But the volume of the latter is a problem we have already solved; it equals $|\text{det} D \varphi(u)| \cdot \text{vol }U$. So what we have "shown" is if $U \subset A$ is a small enough rectangle, then \begin{equation} \text{vol }\varphi(U) \approx \text{vol }D \varphi(u)[U] = |\text{det } D \varphi(u)| \cdot \text{vol }U \end{equation}
So if we now "sum over" all such small open rectangles $U$ which cover $A$, it should make sense why the equation
\begin{equation} \int_{\varphi(A)} 1 \, dv = \int_A |\text{det}\, D \varphi(u)| \, du \end{equation} is true.
Conclusion: The key geometric insight is that for a rectangle $U$ and a linear transformation $\varphi$, the equation \begin{equation} \text{vol }\varphi(U) = |\text{det }\varphi| \cdot \text{vol }U \end{equation} is true. The rest of the proof of change of variables comes by approximating a non-linear $\varphi$ linearly using its derivative $D \varphi(u)$. This is of course nowhere near rigorous, but making these approximations and vague statements more precise is exactly the purpose of the change of variables theorem.
Let $v,\,\varphi$ be $n$-dimensional. We start from a multivariate chain rule $dv_i=\sum_jJ_{ij}d\varphi_j,\,J_{ij}:=\frac{\partial v_i}{\partial\varphi_j}$. Some multilinear scalar function $f(J)$ of the $n\times n$ matrix $J$ satisfies $d^nv=f(J)d^n\varphi$.
Well, the determinant is very special.
In the special case where $v_i=\lambda_i(\varphi_i)$, univariate integration by substitution gives $f(J)=\prod_i|\lambda_i^\prime|=|\det J|$ with $J=\operatorname{diag}(\lambda_i)$. Since we can apply successive variable changes, $f(J_1J_2)=f(J_1)f(J_2)$. The only $f$ that'll do the job is $f(J):=|\det J|$.
(To finish a proof of this, you could use differential forms to show $f$ is antisymmetric in exchanging rows of $J$, as is $\det J$, because $dv_idv_j=-dv_jdv_i$. But we prefer to have $dv_idv_j=dv_jdv_i$ viz. Fubini's theorem, hence the modulus signs.)
