1D Problem and LF Solver

We have 
	({{Q, =, (({{h}, {h u}, {h v}}))}})
The PDE is 
	({{(Subscript[Q, t])+(Subscript[(F(Q)), x])+(Subscript[(G(Q)), y]), =, H(Q)}})
where 
	({{F(Q), =, (({{(Subscript[Q, 2])}, {(((Subscript[Q, 2])^2)/(Subscript[Q, 1]))+(1/2)g ((Subscript[Q, 1])^2)}, {(((Subscript[Q, 2 ])(Subscript[Q, 3]))/(Subscript[Q, 1]))}}))}}) and ({{G(Q), =, (({{(Subscript[Q, 3])}, {(((Subscript[Q, 2 ])(Subscript[Q, 3]))/(Subscript[Q, 1]))}, {(((Subscript[Q, 3])^2)/(Subscript[Q, 1]))+(1/2)g ((Subscript[Q, 1])^2)}}))}}) 
with 
	({{H(Q), =, (({{0}, {-g (Subscript[Q, 1])(Subscript[B, x])}, {-g (Subscript[Q, 2])(Subscript[B, y])}}))}}).

To make stuff simpler I am going to make this a 1D spatial dimension problem.  If nothing depends on y and the y velocity v=0 then we get 
	({{q, =, (({{h}, {h u}}))}})
The PDE is 
	({{(Subscript[q, t])+(Subscript[(f(q)), x]), =, h(q)}})
where 
	({{f(q), =, (({{(Subscript[q, 2])}, {(((Subscript[q, 2])^2)/(Subscript[q, 1]))+(1/2)g ((Subscript[q, 1])^2)}}))}}) 
with 
	({{h(q), =, (({{0}, {-g (Subscript[q, 1])(Subscript[B, x])}}))}}).

Information Velocities

The PDE is 
	({{(Subscript[q, t])+f'(q)(Subscript[q, x]), =, h(q)}})
where 
	({{f'(q), =, (({{0, 1}, {-(((((Subscript[q, 2])/(Subscript[q, 1]))))^2)+g (Subscript[q, 1]), (Subscript[q, 1])}}))}}) .
There are really two "information" speeds encoded in this matrix. They are the eigenvalues!

Df={
{0,1},
{-((((q2/q1)))^2)+g q1,q1}};
λ=Simplify[Eigenvalues[Df]]

Why are these information speeds?
Suppose the matrix was diagonal! Then we would have two decoupled PDES and the speeds would be the diagonal entries.

Clear[g,q1,q2]
{g,q1,q2}={9.8, 1.2,0.1};
{λ,R}=Eigensystem[Df]

The associated eigenvectors give the information that is propagated at that speed.

LF Scheme

Setting up some functions for the flux, bottom, and initial conditions.

Clear[g,a,b,f,q1,q2,B,h,h0,u0]
g=9.8;
{a,b}={0.0,2π};
f[{q1_,q2_}]:={q2,((q2^2)/q1)+0.5 g (q1^2)}
B[x_]:= 1+ 0.1Sin[x];
h[{q1_,q2_},x_]={0, -g q1 B'[x]};
h0[x_]= 2-0.1B[x];
u0[x_]=0;
Plot[{B[x],B[x]+h0[x]},{x,a,b},
PlotRange->{0,4.5}]

Setting up the initial conditions with a given spatial discretization.

Clear[n,Δx,i, q0]
n=42;
Δx=(((b-a))/n);
x[i_]:= a + i Δx;
q0 = Table[ {h0[x],h0[x] u0[x]},{x,a+Δx,b,Δx}];
ListPlot[q0[[All,1]],
PlotLabel->"h at t=0",
DataRange->{a+Δx,b},
Frame->True,
Axes->False,
FrameLabel->{"x","h"}
]

Picking a small Δt (that we hope satisfies the CFL condition) and running for NumSteps steps using the LF algorithm in the paper.  First one step!

Clear[Δt,qOld, qNew]
Δt=0.03;
NumSteps = 800;
qOld=q0;
qNew=q0;
Do[
qNew[[1]]=0.5*(qOld[[n]]+qOld[[2]]) - ((Δt)/(2 Δx))(f[qOld[[2]]]-f[qOld[[n]]])+Δt *h[qOld[[1]],x[1]];
qNew[[n]]=0.5*(qOld[[n-1]]+qOld[[1]]) - ((Δt)/(2 Δx))(f[qOld[[1]]]-f[qOld[[n-1]]])+Δt *h[qOld[[n]],x[n]];
Do[
qNew[[i]]=0.5*(qOld[[i-1]]+qOld[[i+1]]) - ((Δt)/(2 Δx))(f[qOld[[i+1]]]-f[qOld[[i-1]]])+Δt *h[qOld[[i]],x[i]],
{i,2,n-1}];
qOld=qNew,
{NumSteps}]

ListPlot[{q0[[All,1]],qNew[[All,1]]},
PlotLabel->StringForm["t=``",NumSteps*Δt],
PlotStyle->{Red, Green}]

Lets think about the CFL condition for a second.  The information speeds (remember there are two of them) are the eigenvalues of the appropriate derivative of f.  They vary a little bit along the vector but not that much.

evals[{q1_,q2_}]=Eigenvalues[ D[f[{q1,q2}],{{q1,q2}}]];
ListPlot[fred=Transpose[ Map[evals,qNew]]]

Max[Abs[fred]]

Looks as though they are something like ±4.2.   So the CFL condition says that we need to have the discretization "speed" ((Δx)/(Δt)) greater than 4.2.  It seems to do a pretty good job of calculating the onset of instability.

((Δx)/(Δt))