Consider the spherical advection problem: describing the conservation of a property $u$ in a closed spherical domain.
$$ \frac{\partial u}{\partial t}+\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2u\ \boldsymbol{v} \right)=0 $$
In the domain $0<r<1,\ t\geq0$, with the following initial/boundary conditions:
$$ \left.\frac{\partial u}{\partial r}\right|_{r=a}=0, \\ u(r,t=0)=u_i $$ Where $\boldsymbol{v}=\boldsymbol{v}(r,t)$ and $u_i$ are known. ($\boldsymbol{v}\leq0$). Other posts (here, and here) suggested discretizing the conservative form. Forward explicit (downwind) discretization leads to: $$ r\to i\Delta r,\ t\to n\Delta t,\ u(r,t)\to u_i^n \\ \frac{\partial u}{\partial t} \to \frac{u^{n+1}_i-u^{n}_i}{\Delta t} \\ \frac{\partial(r^2u\ \boldsymbol{v})}{\partial r} \to \frac{ ((i+1)\Delta r)^2 u_{i+1}^n \boldsymbol{v}_{i+1}^n- (i\Delta r)^2 u_{i}^n \boldsymbol{v}_{i}^n }{\Delta r} $$ And the equation becomes: $$ \frac{u^{n+1}_i-u^{n}_i}{\Delta t}= -\frac{1}{(i\Delta r)^2} \frac{ ((i+1)\Delta r)^2 u_{i+1}^n \boldsymbol{v}_{i+1}^n- (i\Delta r)^2 u_{i}^n \boldsymbol{v}_{i}^n }{\Delta r} $$
My problems are:
- I cannot use this equation at $r=1$ (last node) because i have no boundary condition to fix the ghost node at $r=1+\Delta r$.
- I'm unsure how to write the equation at $r=0$. I can either use L'Hôpital to write $$ \frac{\partial u}{\partial t}+ \lim_{r\to0}\left(\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2u\ \boldsymbol{v} \right)\right) = \frac{\partial u}{\partial t}+ \lim_{r\to0}\left(\frac{1}{2r}\frac{\partial^2 }{\partial r^2}\left(r^2u\ \boldsymbol{v} \right)\right) =\\ \frac{\partial u}{\partial t}+ \frac{1}{2}\frac{\partial^3 }{\partial r^3}\left(r^2u\ \boldsymbol{v} \right) =0 $$ And discretize the above or I can use the forward derivative and write $$ \frac{u_1^{n}-u_0^{n}}{\Delta r}=0\therefore u_0^{n}=u_1^{n} $$
What is the best way to discretize the problem and how do I address the first and last nodes?
Edit:
I found another post (here) with almost the same problem. One of the answers suggested an interpolation for the stencil at the last node, which i tried with unsuccessful results. The solution returns values greater than $1$ at $r=1$ and does not converge at all as the nodes increase (unstable). The image below shows the solutions for 21, 51 and 101 spatial nodes. The CFL number is below $1$. Solution should stabilize at $0.91$, and be less than $1$ for every instant.
