What does divergence of scalar times vector vector field physically mean? We know that:
$\nabla \cdot (f \vec{A}) = f \nabla \cdot \vec{A} + \vec{A}\cdot(\nabla f)$
Now divergence of any vector field can be understood in terms of whether the concerning flux is outgoing ($\nabla \cdot \vec{A} < 0$) or incoming ($\nabla \cdot \vec{A} > 0$). If your $\vec{A}$ is velocity field, then its divergence represents the change in volume.
From above equation, we can see that $\nabla \cdot (f \vec{A})$ depends upon (sign) of scalar field: $f$ and also its gradient. Can someone help me to understand how we can physically interpret the above equation?
My understanding is that if $\nabla \cdot(f \vec{A})$ < 0 in region $R$, then it implies scalar field $f$ is getting removed from region $R$, and if  $\nabla \cdot (f \vec{A})$ > 0, then scalar field $f$ is adding to region $R$. I reached to this conclusion by thinking of a continuity equation and taking $f = \rho$ (density) and $\vec{A} = \vec{v}$ (velocity field)
 A: The physical meaning of $\nabla (f \vec A)$ is the same as for a single vector field, namely it is a measure of a flow out of a region. What exactly it is that flows depends on what quantities are described by $f$ and $\vec A$. Taking your example of $f=\rho$ being a mass density and $\vec A = \vec v$ being a velocity field, $f \vec A = \rho \vec v$ is just the mass current density, which I will call $\vec j$. Then
$$
\nabla (f \vec A) = \nabla(\rho \vec v) = \nabla \vec j~,
$$
and the physical interpretation of the last expression should be clear.
Now let's look at your expansion
$$
\nabla (f \vec A) = (\nabla f) \vec A + f \nabla \vec A~.
$$

*

*If $f$ is constant, the first term vanishes and we get $\nabla(f \vec A) = f \nabla \vec A$, in agreement with $\nabla$ being $\mathbb C$-linear. The physical interpretation of this is that both the vector field and its divergence get scaled by some constant $f$.

*If $f$ is not constant, but differentiable, it can be approximated around any point p as
$$
f(x) = f(p) + \underbrace{f'(p) (x-p)}_{=(\nabla f(p)) (x-p)} + \mathcal O((x-p)^2)~.
$$
The constant term $f(p)$ leads to the appearance of $f \nabla \vec A$ in $\nabla (f \vec A)$ also for non-constant $f$. The physical interpretation of this term is the same as for $f$ being constant, only that there is potentially a different constant for every point in space.

*If $f$ is not constant, there is, however, another possible reason why $f\vec A$ would change: supposing $f$ increases in the direction of $\vec A$, then the absolute value of the flow $|f \vec A|$ increases and there is more flow out of the region in question than into it (positive divergence). In this case, $\nabla f$ has a component parallel to $\vec A$, so $(\nabla f) \vec A > 0$. This motivates that $(\nabla f) \vec A$ indeed contributes to $\nabla (f \vec A)$.

*Because we only work to linear order when taking derivatives (derivation is linearisation), no interferences between the terms explained in the last two points have to be considered, and they can just be summed up, which yields the equation you asked about.

