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.

• In case anyone wants to look further into the origin of this problem, this is the article I am trying to solve numerically. – eg.Zeta Aug 24 '16 at 11:28
• I'm confused. With $\nu<0$, you have material flowing inward from the outside boundary, so your discretization is upwind-biased (not downwind), and you need a boundary condition at the outside of the sphere, which you haven't provided. Furthermore, it seems intuitively like the exact solution should develop a singularity at the origin, since all the mass arrives there eventually. – David Ketcheson Aug 25 '16 at 10:58