I figured it out. Understand that the substantial derivative exists for only field variables i.e. functions that vary with position vector x and time t.
Thus,
$$
\frac{DM}{DT}
$$
doesn't really make sense.
For a fixed control volume the mass can only vary with time.
Due to the complexity of the topic of the question, it is not possible to expound here all the associated analytic developments, therefore I try to answer by giving precise references and surveying the results presented there.
- Do solutions for $f$ exist?
Yes: to see this, note that, after formally manipulating your equation (perhaps going backward in the steps of your deduction), we get
$$
\begin{split}
\nabla\cdot\nabla f(\vec{r})&=\left( \frac{1}{2}\vec{B}\times\vec{r} - \nabla f(\vec{r}) \right)\cdot\frac{\nabla \rho(\vec{r})}{\rho(\vec{r})}\\
\rho(\vec{r})\nabla\cdot\nabla f(\vec{r})&=\left( \frac{1}{2}\vec{B}\times\vec{r} - \nabla f(\vec{r}) \right)\cdot\nabla \rho(\vec{r})\\
\rho(\vec{r})\nabla\cdot\nabla f(\vec{r})&=\left( \frac{1}{2}\vec{B}\times\vec{r} - \nabla f(\vec{r}) \right)\cdot\nabla \rho(\vec{r})\\
\rho(\vec{r})\nabla\cdot\nabla f(\vec{r})+\nabla \rho(\vec{r})\cdot\nabla f(\vec{r})&=\frac{1}{2}\vec{B}\times\vec{r}\cdot\nabla \rho(\vec{r})\\
\end{split}
$$
i.e.
$$
\nabla\cdot\big(\rho(\vec{r})\nabla f(\vec{r})\big)=\frac{1}{2}\vec{B}\times\vec{r}\cdot\nabla \rho(\vec{r}).\label{1}\tag{1}
$$
Equation \eqref{1} is a standard linear elliptic equations in divergence form, and as such it admits a fundamental solution (and also a Green function) i.e. a solution of the equation
$$
\nabla\cdot\big(\rho(\vec{r})\nabla \mathscr{E}(\vec{r},\vec{s})\big)=\delta(\vec{r}-\vec{s})\label{2}\tag{2}
$$
where $\delta$ is the usual Dirac distribution, under fairly general conditions. In particular $\rho(\vec{r})$ is required only to be a bounded and measurable function: moreover, the problem has a solution even if, instead of a scalar function $\rho(\vec{r})$, we have a second-order tensor function
$$
\vec{r}\mapsto\overline{\overline{\boldsymbol{\rho}}}(\vec{r})\in\Bbb R^{3\times 3}\quad\vec{r}\in\Bbb R^3 \label{3}\tag{3}
$$ with bounded measurable components. This result is due to Walter Littman, Guido Stampacchia and Hans Weinberger: the paper [2] (where the original result is to be found) is not an easy read and even the used notation is not a modern one, while the set of lecture notes [3] is a better readable source (with also an updated notation) but is written in French. A more modern reference is [1], where the authors prove also that the tensor of the coefficients \eqref{3} can be asymmetric.
- Is it possible to say something else than the above equation on the form of the solutions $f$? (Especially: is there an analytical solution?)
Yes (with some cautions): for the general case, several properties of the solution to equation \eqref{2}, and these properties "interact" with the ones of the given function $\frac{1}{2}\vec{B}\times\vec{r}\cdot\nabla\rho(\vec{r})$ (have again a look at the items in the bibliography for the details and also at this answer). For example it is known that the quotient between $\mathscr{E}(\vec{r},\vec{s})$ and the fundamental solution of Laplace equation is bounded by two positive constants $K$ and $K^{-1}$. And also you can use $\mathscr{E}(\vec{r},\vec{s})$ to construct a fairly explicit general solution of \eqref{1} as a convolution integral:
$$
f(\vec{r})=\int\limits_{\Bbb R^3}\mathscr{E}(\vec{r},\vec{s})\big(\vec{B}\times\vec{s}\cdot\nabla\rho(\vec{s})\big)\mathrm{d}\vec{s}
$$
However, it is not always possible to construct explicitly the fundamental solution $\mathscr{E}(\vec{r},\vec{s})$ (in this respect see also this answer) thus, apart from the case $\rho(\vec{r})\equiv\mathrm{const}.$ where \eqref{1} reduce to the Laplace equation
$$
\nabla\cdot\nabla f(\vec{r})=0,
$$
I am not aware of the existence of any exact (meant as expressed by means of more or less elementary functions) solutions to \eqref{1}.
Bibliography
[1] Michael Grüter and Kjell-Ove Widman (1982), "The Green function for uniformly elliptic equations", Manuscripta Mathematica 37, pp. 303-342 Zbl 0485.35031.
[2] Walter Littman, Hans Weinberger and Guido Stampacchia (1962), "Regular points for elliptic equations with discontinuous coefficients", Annali della Scuola Normale Superiore di Pisa, Classe di Scienze, serie III, Vol. 17, n° 1-2, pp. 43-77, MR161019, Zbl0116.30302.
[3] Guido Stampacchia (1966), "Équations elliptiques du second ordre à coefficients discontinus" (notes du cours donné à la 4me session du Séminaire de mathématiques supérieures de l'Université de Montréal, tenue l'été 1965), (in French), Séminaire de mathématiques supérieures 16, Montréal: Les Presses de l'Université de Montréal, pp. 326, ISBN 0-8405-0052-1, MR0251373, Zbl 0151.15501.
Best Answer
This is basically a variable coefficient transport equation. The method of characteristics is well-suited for such equations.
Altering the notation a bit, you have $$\partial_tu + (x\partial_x - y\partial_y)u = 0.$$ The idea is to transform the above PDE into an ODE along suitable curves. Taking $u=u(x(s),y(s),t(s))$, we see that we must have $$\frac{dx}{ds} = x(s), \ \frac{dy}{ds}=-y(s), \ \frac{dt}{ds}=1, \ \frac{du}{ds}=0.$$ Assuming $t(0)=0$, we get $t=s$. Further, letting $x(0)=x_0$ and $y(0)=y_0$, we get $$x=x_0e^{s} = x_0e^{t} \text{ and } y=y_0e^{-s} = y_0e^{-t}.$$ By letting $u(0) = \varphi(x_0,y_0)$, we get $u(x,y,t) = \varphi(xe^{-t},ye^{t},t)$. However, we are given $u(x,y,0) = x^2 + y^2$. This implies that $$u(x,y,t) = x^2e^{-2t}+y^2e^{2t}.$$