Fix shallow water equation

[xtalax]

Instability found in following solution, what could be causing this?

using ModelingToolkit, MethodOfLines, DomainSets, OrdinaryDiffEq
using Plots
using IfElse

# Define the system of equations
@parameters t x
@variables u(..) η(..)
#hu(..)

Dt = Differential(t)
Dx = Differential(x)

tmin = 0.0;
tmax = 6000.;

xmin = 0.0;
xmax = 200.;

dx = 2.0;

slope = 0.001;
n= 0.03;
order= 2;


elevation(x) =  (xmax - x) * slope;
h(t,x) = η(t,x) - elevation(x);


source(t) = ifelse(t < 3600, 1/60/1000, 0.0)
@register_symbolic source(t)
g = 9.81;

exp = -4/3

eqs = [
    Dt(η(t,x))  ~ source(t) - u(t,x) * Dx(η(t,x)) + u(t,x) * (-slope) - h(t,x) * Dx(u(t,x)) ,
    Dt(u(t,x))  ~ - g * Dx(η(t,x))  - u(t,x) * Dx(u(t,x))  - g * u(t,x)^2 * n^2 * IfElse.ifelse(h(t,x) > 1e-5, h(t,x)^exp, (1e-5)^exp)
]

domain = [x ∈ Interval(xmin,xmax),
        t ∈ Interval(tmin,tmax)]


bcs = [
    η(tmin,x)     ~ elevation(x)+1e-5,
    u(tmin,x)     ~ 0.0,
    η(t,xmin)     ~ elevation(xmin),
    u(t, xmin)    ~ 0.0,
]

@named pdesys = PDESystem(eqs, bcs, domain, [t, x], [u(t,x),  η(t,x)])
discretization = MOLFiniteDifference([x => dx], t, jac = true, sparse = true) #
grid = get_discrete(pdesys, discretization)

prob = ModelingToolkit.discretize(pdesys, discretization)
@time sol = solve(prob, Tsit5())

anim = @animate for (i, T) in enumerate(sol.t)
    plot(map(d->sol[d][i], grid[u(t,x)]), label="u", title = "t = $T")
    plot!(map(d->sol[d][i], grid[η(t,x)]), label="η")
end
gif(anim, "plots/shallow_water.gif", fps = 10)

Plot:

[YichengDWu]

Looks oscillatory to me. WENO methods may fix this.

[xtalax]

Fixed with WENO and implicit methods