Spherical laplacian example seems to give moderate error between exact and numerical solution

[00krishna]

Hello. I was testing out some code from the MOL test cases for a spherical laplacian problem--basically a problem using polar coordinates. Here is a link to the test case simulation seems to give a relatively high level of error, so I was not sure what the problem was. I tried to increase the order of discretization and some other tricks, but that did not solve the issue. Below is a complete example and a copy of a GIF file that shows the problem.

I am not sure if this level of error is acceptible, or more than acceptible? The error does not seem to be unstable or growing, but there is a clear difference between the analytic and numerical solutions.

using ModelingToolkit, MethodOfLines, LinearAlgebra, OrdinaryDiffEq, DomainSets
using ModelingToolkit: Differential


#=
Basic problem: Spherical Laplacian

In this problem, the Laplacian is given in polar coordinates. So we can 
look at diffusion in polar coordinates. 

Du_t = 1/r^2 *  d/dr[(r^2 * Du_r)]

This is again a manufactured solution.
=#

u_exact(r, t) = exp.(-t) * sin.(r) ./ r

@parameters t r
@variables u(..)
Dt = Differential(t)
Dr = Differential(r)

# equation system 

eq = Dt(u(t, r)) ~ 1/r^2 * Dr(r^2 * Dr(u(t, r)))

bcs = [
        u(0, r) ~ sin(r)/r,
        Dr(u(t, 0)) ~ 0,
        u(t, 2) ~ exp(-t) * sin(1)
]

domains = [
            t ∈ Interval(0.0, 2.0),
            r ∈ Interval(0.0, 2.0)
]

@named pdesys = PDESystem(eq, bcs, domains, [t, r], [u(t,r)])

# discretize

dr = 0.1
order = 4
discretization = MOLFiniteDifference([r => dr], t, approx_order=order)

prob = discretize(pdesys, discretization)

sol = solve(prob, Rodas4(), saveat=0.1)

# Plot results and compare with exact solution
discrete_r = sol[r]
discrete_t = sol[t]
solu = sol[u(t, r)]

using Plots
plt = plot()

anim = @animate for (i, T) in enumerate(discrete_t)
    exact = u_exact(r, T)
    print(exact)
    plot(discrete_r, exact, seriestype = :scatter, label="analytic")
    plot!(discrete_r, solu[i,:], label="numeric", ylims=(-0.1, 1.1))
    #print(abs2.(exact-solu[i,:]))
end

gif(anim, "mol_diffusion_1d_spherical_07.gif", fps=5)

Finally here is a copy of the plot GIF.