Analytical treatment
Using spherical polar coordinates
$$
x = r \cos(\phi) \sin(\theta)
\quad ; \quad
y = r \sin(\phi) \sin(\theta)
\quad ; \quad
z = r \cos(\theta)
$$
one obtains the most convenient analytical formulation of the problem.
In all subsequent cases, the drift velocity $\vec{v}_{dr}$ is zero and the
velocity of the solar wind is radial & constant: $\vec{v}_{sw}=v_0 \vec{e}_r$.
Hence the term $-\gamma (\nabla\cdot \vec{v}_{sw})$ can be evaluated as
$-\gamma v_0 \nabla\cdot \vec{e}_r = -\gamma v_0 (2 / r) $.
Here comes a compilation of (fictive) one-dimensional analytical solutions
of the CR equation. It is remarked that this version of the theory
actually matches reasonably well with (parts of) the numerical solutions.
So I have confidence that both theory and code are correct, at last.
MODE=1: $\; \kappa = radial = r \; ; \; v_0 = 0 $
The Partial Differential Equation (P.D.E.) is reduced, for $\kappa$ equal to
the scalar $r$, and for $v_0$ equal to zero, to the following equation:
$$ \frac{d}{dr} r \frac{dp}{dr}+2\frac{dp}{dr} = 0 $$
$$ \left( \frac{d}{dr}+\frac{2}{r} \right) r \frac{d}{dr} p = 0 $$
Using a basic formula of Operator Calculus
$\; \large \frac{d}{dx}+f(x) = e^{-\int f(x)dx }\,\frac{d}{dx}\,e^{+\int f(x)dx } \;$
with $ \; \exp( \int 2/r \, dr ) = r^2 $ :
$$ \left( \frac{1}{r^2} \frac{d}{dr} r^2 r \right) \frac{d}{dr} p = 0 $$
The general solution is $ p = C / r^2 $. Boundary
conditions are chosen spherically symmetric: $p(RI)=0$ and $p(RO)=1$.
Then the solution is: $$ p(r)=\frac{1/RI^2-1/r^2}{1/RI^2-1/RO^2} $$ It has been
seen that this solution agrees VERY well with the numerical values.
MODE=2: $\; \kappa = unity \; ; \; v_0 = 0 $
The P.D.E. is reduced, for $\kappa$ equal to the unity tensor, and for $v_0$
equal to zero, to the Laplace equation in spherical coordinates. The boundary
conditions are also chosen spherically symmetric: $p(RI)=0$ and $p(RO)=1$.
Then the solution is: $$ p(r)=\frac{1/RI-1/r}{1/RI-1/RO} $$ Since it is well
known that the potential field of a point source goes like $1/r$. It has been
seen that this solution also agrees VERY well with the numerical values.
MODE=3: $\; \kappa = radial = r $
The P.D.E. is reduced to an ordinary differential equation: $$
\frac{d}{dr}r\frac{dp}{dr}+(2-v_0)\frac{dp}{dr}-\gamma v_0\frac{2}{r}p = 0 $$
Here $v_0=0.5984$ and $\gamma=1.663151$.
This is a differential equation of type Euler in $p(r)$, $a$ and $b$ constant:
$$ r^2 p" + (a+1) r p' + b p = 0 $$
Where $ a=2-v_0 $ ; $ b=-2 \gamma v_0 $.
The general solution is, in case $\lambda_1 \neq \lambda_2 $ :
$$ p(r) = C_1 r^{\lambda_1} + C_2 r^{\lambda_2} \qquad
C_1,C_2 \mbox{ arbitrary } $$
$\lambda_{1,2}$ are the solutions of the characteristic equation:
$$ \lambda^2 + a \lambda + b = 0 $$
Substituting herein the numerical values of $v_0$ and $\gamma$ yields:
$$ \lambda_1 = 0.8745031 \qquad ; \qquad \lambda_2 = -2.2761030 $$
The constants $C_{1,2}$ are determined by the boundary conditions, as before.
The result is:
$$ p(r)=\frac{r^{\lambda_1}RI^{\lambda_2}- r^{\lambda_2}RI^{\lambda_1}}
{RO^{\lambda_1}RI^{\lambda_2}-RO^{\lambda_2}RI^{\lambda_1}} $$
Which also matches VERY well with the numerical solution.
MODE=4: $\; \kappa = unity = 1 \; ; \; \gamma = 0 $
If we assume that $v_0 \neq 0$, and specialize for $\gamma = 1$, then
CR reduces to the following ordinary differential equation: $$
\frac{d^2p}{dr^2}+(2/r-v_0)\frac{dp}{dr}- v_0\frac{2}{r}p = 0 $$
The solution of this differential equation can be found again by employing
Operator Calculus
(which is kind of routine for me if I want to find "exact" solutions):
$$ \left( \frac{d}{dr} + \frac{2}{r} \right)
\left( \frac{d}{dr} - v_0 \right) p = 0 $$ Using the basic formula
$ \; \large \frac{d}{dx}+f(x) = e^{-\int f(x)dx }\,\frac{d}{dx}\,e^{+\int f(x)dx } \;$ :
$$ \left( r . \frac{1}{r^2} \frac{d}{dr} r^2 \right)
\left( e^{v_0 r} \frac{d}{dr} e^{-v_0 r} \right) p = 0 $$
Systematic integration gives:
$$ \frac{d}{dr}e^{-v_0 r} p = C_1 \frac{e^{-v_0 r}}{r^2} $$
$$ p(r)= C_1 e^{v_0 r} \int \frac{ e^{-v_0 r}}{r^2} dr + C_2 e^{v_0 r} $$
So there exist the following two elementary solutions:
$$ p_1(r)=e^{v_0 r} \int \frac{ e^{-v_0 r}}{r^2} dr $$ $$ p_2(r)= e^{v_0 r} $$
The integral can be worked out further by partial integration, with $t=-v_0r$:
$$ \int \frac{e^t}{t^2}dt = - \frac{e^t}{t} + \int \frac{e^t}{t}dt $$
Giving for the first solution:
$$ p_1(r) = \frac{1}{v_0 r} + e^{v_0 r}.Ei(- v_0 r) $$
$$ Ei(x)= \int_{-\infty}^{x} \frac{e^t}{t} dt $$
If you don't believe this, you can substitute back and check out with MAPLE:
p:=1/(v0*r)+exp(v0*r)*Ei(-v0*r);
r*diff(diff(p,r),r)+(2-v0*r)*diff(p,r)-2*v0*p;
simplify(");
quit;
Together with the boundary conditions, the solution finally becomes:
$$ p(r)=\frac{p_1(RI)p_2(r)-p_2(RI)p_1(r)}{p_1(RI)p_2(RO)-p_2(RI)p_1(RO)} $$
Which in turn can be compared with numerical results. The latter is somewhat
difficult for two reasons: the exponential integral function Ei(x) is unknown
to Fortran, and the analytical result has a very steep gradient for the nominal
value of $v_0$. The first problem can be adressed by invoking the NAG Fortran
Library: the exponential integral function Ei(x) is calculated by a function
S13AAF(X,IFAIL). The second problem can be adressed by lowering the value of
$v_0$ by a factor of 100.
MODE=5: $\; \kappa = unity = 1 \; ; \; \gamma = 0 $
One more analytical solution of the CR equation can be constructed, resulting
in a total of 5 exact solutions so far.
If we assume that $\kappa$ is the unity tensor, $v_0 \neq 0$, and specialize
for $\gamma = 0$ instead of $\gamma = 1$, then the CR Confusion (:-) P.D.E.
reduces to the following ordinary differential equation: $$
\frac{d^2p}{dr^2}+(2/r-v_0)\frac{dp}{dr} = 0 $$
The solution of this differential equation can be found again by
Operator Calculus :
$$ \left[ \frac{d}{dr} + \left( \frac{2}{r} - v_0 \right) \right]
\frac{d}{dr} p = 0 $$ Using the basic formula
$ \frac{d}{dx}+f(x) = e^{-\int f(x)dx }\,\frac{d}{dx}\,e^{+\int f(x)dx } \,$ :
$$ \left( \frac{1}{r^2} e^{v_0 r} \frac{d}{dr} r^2 e^{- v_0 r} \right)
\frac{d}{dr} p = 0 $$
Systematic integration gives:
$$ p(r)= C_1 \int \frac{ e^{v_0 r}}{r^2} dr + C_2 $$
The integral can be worked out as in the previous case, giving the solution,
apart from constants:
$$ p_1(r) = - \frac{e^{v_0 r}}{v_0 r} + Ei(v_0 r) $$
If you don't believe this, you can substitute back and check out with MAPLE:
p:=-exp(v0*r)/(v0*r)+Ei(v0*r);
diff(diff(p,r),r)+(2/r-v0)*diff(p,r);
simplify(");
quit;
Together with the boundary conditions, the solution finally becomes:
$$ p(r)=\frac{ p_1(RI)-p_1(r) }{ p_1(RI)-p_1(RO) } $$
Which in turn can be compared with numerical results. It should be emphasized,
however, that the analytical, and hence also the numerical solution, is
only known for the case $v_0 < 0$. It is just zero in other cases!
Implementation
The five analytical solutions have been implemented for testing purposes as
the Fortran function EXAKT(R), where R is a radial distance from the origin.
The solution is selected by an index MODE, which is ranging from 1,..,5 and
resides in the common block /TEST/. MODE is also the index which determines
what coefficients to calculate in subroutine COEFFS. Needless to say that
the choice of the exact solution and the choice of the coefficients should
lead to a match between the exact and the numerical solutions respectively.
Under normal circumstances, COEFFS contains the entire physics of the problem.
Given the convection-difusion equation as such, together with its boundary
conditions and the geometry of the calculation domain, only the coefficients
are yet to be determined, by proper physical modelling. The relationship
between the parameters of COEFFS and the (test) physics of CR is as follows:
$$
\stackrel{\leftrightarrow}{\kappa} = \left[ \begin{array}{ccc}
a_{xx} & a_{xy} & a_{zx} \\ a_{xy} & a_{yy} & a_{yz} \\ a_{zx} & a_{yz} & a_{zz}
\end{array} \right] \qquad = \mbox{symmetric}
$$ $$
- v_0 \frac{x}{r} = a_x \qquad - v_0 \frac{y}{r} = a_y \qquad
- v_0 \frac{z}{r} = a_z
$$ $$
- \gamma v_0 \frac{2}{r} = a_0
$$
With our test examples, the $\kappa$ tensor can be unity or radial.
In the latter case, the coefficients are calculated by:
$$
a_{xx} = \frac{x^2}{r} \qquad a_{yy} = \frac{y^2}{r}
\qquad a_{zz} = \frac{z^2}{r} $$
$$ a_{xy} = \frac{xy}{r} \qquad a_{yz} = \frac{yz}{r}
\qquad a_{zx} = \frac{zx}{r} $$